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1.  Satellite  Orbits 

> 

' 

Initiator:  E.  Robinson 

Project  No:  0001  Problem  No:  4517 

I ' 

| • A broad  range  of  research  analyses  conducted  at  AFGL  requires,  as  in** 

\ 

put,  precise  satellite  orbit  data.  A researcher  might,  for  example, 

be  interested  in  using  on  board  sensors  to  map  a particular  physical 

j 

quantity  through  space.  He  would  need  to  know  the  satellite's  position 
as  a function  of  time  in  order  to  interpret  telemetry  data,  or  he  might 
need  to  know  the  satellite's  precise  location  relative  to  some  spatial- 
ly distributed  physical  quantity  such  as  solar  radiation  at  the  edge  of 
the  earth's  shadow.  Programs  LOKANGL,  ROPP  as  well  as  other  adapta- 

i 

tions  of  conventional  orbit  determination  programs  are  available  to 
provide  the  needed  data.  Each  has  Its  own  unique  features  in  terms  of 
the  type  of  input  data  on  which  calculations  are  based,  and  methods  of 
treating  perturbations  to  orbital  dynamics  such  as  drag  and  the  higher 
order  terms  of  the  geopotential.  Efforts  are  continuously  being  con- 
ducted to  evaluate  and  Improve  the  accuracy  of  these  programs  and  to 
generalize  and  extend  their  applicability.  These  studies  are  reported 
i in  the  following  subsections. 
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1.1  LOKANGL  Features 

\ 

Program  LOKANGL  has  continued  to  be  used  and  developed  for  most  AFGL  > 

satellite  ephemerls  generation  and  study  purposes.  Improvements  have 
been  made  in  providing  satellite  observation  data  for  flying  aircraft, 
and  in  the  satellite  solar-illumination  routine.  Additional  capabi- 
lities were  Implemented  in  order  to  handle  satellite  thrusting;  the  ’ 

perigee-apogee  determination  algorithm  was  also  modified  to  function 
properly  under  thrust  adjustments.  Orbital  elements  are  now  accepted 
in  any  of  five  formats:  a)  SCF  position- velocity  sets  In  subroutine 
DELEM;  b)  2 -card  ADC,  3-card  ADC,  and  osculating  Keplerian  sets  in 
subroutine  CONVRT;  c)  original  5-card  ADC  Keplerian  sets  directly  In 
LOKANGL. 

I 

Output  options  include  print-out  of  selectable  station-independent  or 
station-related  parameters;  and  a binary  ephemerls  tape  that  is  in  the 
original  format  defined  in  the  ancestral  source  program1  and  that  is 
retained  with  minor  changes.  The  current  Input  format  is  presented  in 
Figure  1;  and  the  binary  tape  format  Is  described  In  Figure  2. 

Program  LOKANGL  Is  designed  to  convert  any  of  the  element  sets  provided 
into  standard  Keplerian  form,  and  then  to  generate  an  ephemerls  after 
revising  element  rates  of  change  where  Interpolation  Is  possible.  The 
effects  of  first  order  geopotential  and  atmospheric  perturbations  only 
are  Included.  Sources  of  Inaccuracy  In  the  trajectory  may  be  due  to 
element  set  discrepancies  Introduced  by  the  orbit  determination  pro- 
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Figure  2.  Binary  TAPE3  Output  Format 
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I 

1 

I 

grams,  or  due  to  the  simplifications  of  the  LOKANGL  model.  A third 

source  of  error  entails  the  Interpretation  of  element  sets  as  provided 

by  various  agencies.  Position- velocity  sets  are  definitive  In  this 

2 

respect,  and  an  earlier  report  described  studies  and  conversions  that 

i 

were  employed  to  standardize  results  when  using  mean  elements.  Special 
studies  conducted  to  evaluate  errors  due  to  the  LOKANGL  simplifications 
are  summarized  In  Section  1.2,  and  show  that  less  than  2 sec.  In- track 
variations  occur  when  using  dally  element  sets.  Discrepancies  between 
element  sets  are  found  to  have  comparable  In- track  variations,  and  the 
normal  procedure  Is  to  discard  erratic  element  sets  after  analytical  or 
graphical  evaluation  of  results  with  and  without  these  sets. 


f 

I 

I 

1.1.1  Observing  Satellite  From  Aircraft 

Initiator:  J.  Buchau  > 

Project  No:  4643  Problem  No:  4965 

i 

f 

The  capability  of  program  LOKANGL  to  produce  the  position  of  a satel- 
lite with  respect  to  an  observation  station  through  the  use  of  sub- 
routine  CORFL  and  special  service  program  FLTRANS  was  consolidated  by 
the  Incorporation  of  subroutine  FLTRANS  into  LOKANGL  and  the  elimin- 
ation of  special  service  program  FLTRANS.  The  present  modification 
now  allows  the  researcher  direct  access  to  LOKANGL  with  the  aircraft 
position-time  cards  in  their  original  format,  with  minimal  adjustments 

to  the  standard  LOKANGL  input  deck.  Aircraft-satellite  look  angles 

. ' 

continue  to  be  simulated  by  a series  of  moving  stations  representing 
the  aircraft  flight  segment.  LOKANGL  output  has  been  modified  to 
clearly  Indicate  the  flight  parameters  and  flight  header  Information 
of  Interest,  and  output  paper  requirements  have  been  reduced  by  80-85% 
in  this  consolidated  presentation.  Figure  1 reflects  the  updated  Input 
format,  and  Figure  3 presents  sample  aircraft-satellite  look  angle  out- 

! 

put  display  in  Its  present  form. 

Functional  Description 

Subroutine  CORFL  Interpolated  the  aircraft  position  at  any  desired  In- 
stant In  the  given  flight  path.  Since  one  execution  of  LOKANGL  can  now 
supply  the  user  with  the  entire  aircraft-satellite  look  angle  history, 
large  printout  gaps  and  redundancy  In  output  display  Information  have 
been  avoided.  In  the  previous  system,  one  execution  of  LOKANGL  calcul- 
ated for  only  one  of  the  array  of  flight  segments.  Aircraft  position 
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Figure  3b.  Aircraft-LOKANGL  Sample  Printout 


and  time  coordinates  are  written  on  TAPE 7 Immediately  after  they  are 
read  by  subroutine  FLTRANS,  and  subsequent  reading  of  TAPE 7 at  the  ap- 
propriate time  eliminates  the  necessity  of  increased  core  storage,  and 
allows  for  flights  with  a virtually  unlimited  number  of  segments. 

1.1.2  Input  of  Elements  for  Thrusted  Satellites 

Program  LOKANGL  was  extended  to  generate  orbits  for  thrusted  satellites. 
The  standard  elements  are  entered  with  single  thrust  cards  Interspersed 
at  the  appropriate  times.  The  program  was  modified  to  discontinue  use 
of  a prior  element  set  when  the  ephemeris  reaches  a thrust  time,  and 
then  to  project  backward  from  a subsequent  element  set.  Diagnostics 
are  included  to  check  for  Improper  card  sequences.  Capability  for  hand- 
ling two  successive  thrusts  with  no  element  set  In  between  is  feasible, 
but  has  not  been  implemented.  The  system  has  been  used  extensively  with 
SCF  element  sets  to  generate  S3-4  and  PACA  orbits.  Section  1.3  presents 
the  results  of  a short  study  to  determine  the  effect  of  various  thrust 
adjusts  on  S3-4  orbital  parameters. 

1.1.3  Apogee-Perigee  Calculations  For  Thrusted  Satellites 

The  discontinuous  nature  of  the  argument  of  perigee,  the  altitude  at  apo- 
gee, and  the  altitude  at  perigee  necessitated  a revision  of  the  apogee- 
perigee  calculation  procedures.  The  epoch  Is  advanced  progressively  to 
Increase  mean  anomaly  by  180°  and  then  Iteratively  optimized  to  the  time 
of  apogee  or  perigee.  As  long  as  a thrust  does  not  cause  an  adjustment 
of  greater  than  90°  In  the  argument  of  perigee,  the  above  procedure  will 
not  be  confounded.  Figure  4 shows  a standard  one-month  apogee-perigee 
plot  for  S3-4. 

1.1.4  ADC  3-Card  Element  Sets  (C0NVRT3) 

Program  LOKANGL  has  been  modified  to  accept  as  Input  ADC  3-card  element 
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data  sets.  Required  Input  parameters  are  Identical  to  those  supplied 
by  the  ADC  SPADATS  2-card  element  sets,  and  scaling  adjustments  are 
made  by  entry  point  C0NVRT3  of  subroutine  CONVRT  to  the  argument  of 
perigee,  right  ascension  and  year  In  order  to  standardize  values  for 
proper  use  by  LOKANGL.  Figure  1 reflects  the  updated  Input  format, 
and  Figure  5 describes  the  three-card  format  now  acceptable  as  Input 
to  LOKANGL. 


1.2  Orbit  Estimation  Accuracy  Studies 


The  most  significant  variation  between  orbits  generated  by  LOKANGL 
or  ROPP  using  successive  element  sets  shows  up  as  the  In- track  error. 
Due  to  the  earth's  rotation  this  timing  error  also  produces  a ground 
track  error  of  up  to  0.5  km  In  longitude  per  second  of  In-track  time, 
which  clearly  shows  up  when  predicting  high  drag  orbits.  Radial  and 
cross-track  errors  are  generally  small.  The  In- track  difference 
between  orbits  at  a given  epoch  may  be  obtained  from  the  Keplerlan 
elements  at  the  epoch,  using  argument  of  perigee  plus  mean  anomaly 
(W  + M)  as  the  angular  measure  of  In- track  position;  or  from  position- 
velocity  vectors  (T* , Y)  as 


[Pr  " * ^R 


where  the  subscripts  Indicate  reference  (R)  or  computed  (C)  quantities. 


Section  1.2.1  describes  a LOKANGL  Interpolated  orbit  study  which  used 
position-velocity  vector  comparisons,  and  section  1.2.2  presents  a 
ROPP  prediction  study  which  used  osculating  Keplerlan  element  compar- 
isons. 

1.2.1  LOKANGL  Interpolated  Orbits 

Three  hypothetical  satellite  orbits,  each  covering  a 4-day  span,  were 
generated  using  CELEST  with  the  STEM  8x8  geopotential  model  and  the 
Jacchla  64  atmospheric  density  model.  The  orbits  Included  high  drag 
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160  km  perigee  altitude,  moderate  drag  240  km  perigee  altitude,  and 
synchronous  altitude  cases.  The  trajectories  could  be  considered  to 
be  "real"  since  the  CELEST  model  included  considerable  fine  structure. 

Position- velocity  vectors  were  generated  on  an  hourly  basis;  daily 
i vectors  with  estimated  period  decays  were  then  selected  as  the  "SCF"- 

type  elements  for  input  to  LOKANGL,  and  the  resulting  interpolated 

orbit  was  compared  against  the  hourly  vectors.  > 

The  original  approach  taken  In  order  to  match  mean  anomaly  at  the 
element  epochs  was  to  adjust  mean  motion  for  low  drag  orbits,  but  to 
adjust  the  derivative  of  mean  motion  for  high  drag  orbits.  Results 
of  these  tests  showed  that,  using  element  sets  separated  by  4 days 
for  high  drag  orbits,  CELEST-LOKANGL  In-track  differences  exceeded 
100  km  or  15  seconds  of  time,  and  that  this  situation  could  be  im- 
proved only  by  adjusting  both  mean  motion  and  Its  derivative.  The 
procedure  now  Implemented  adjusts  the  derivative  of  mean  motion  based 
on  the  mean  motions  before  adjusting  mean  motion  to  match  mean  anom- 
aly. Maximum  In- track  timing  error  has  thus  been  reduced  to  2-3 
seconds  over  the  4-day  span.  The  low  drag  and  synchronous  orbit 
cases  are  essentially  unaffected  by  this  modification  and  have  main- 
tained less  than  2 seconds  of  error  over  the  4-day  span. 

In  addition  to  the  above  limitations,  elements  transmitted  by  agen- 
cies are  subject  to  orbit  determination  Inaccuracies.  This  problem 
Is  evidenced  by  erratic  elements  for  the  same  epoch  when  transmitted 
on  different  days.  10  km  or  greater  variations  have  been  observed, 
and  therefore  the  absolute  accuracy  of  LOKANGL  ephemerldes  Is  limited 
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to  about  25  km. 
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1.2.2  ROPP  Prediction  Study 

S3-4  orbit  predictions  were  made  over  a 10  day  span  in  April  1978 
using  program  ROPP.  A two-day  thrust-free  span  was  first  employed 
to  determine  an  optimized  drag  factor,  which  was  established  as 
CDAM  = 0.0104796  ft^/lb  for  all  the  runs.  SCF  vector  and  ADC 
2-card  Keplerian  elements  sets  were  used,  and  the  in-track  position, 
as  measured  by  osculating  argument  of  perigee  plus  mean  anomaly,  is 
presented  in  Table  1 (SCF)  and  Table  2 (ADC)  at  integral  hours  close 
after  the  SCF  element  set  epochs.  Various  conclusions  are  summarized 
below. 

1)  Small  initial  differences  in  the  semi -major  axis  cause  the 
in-track  error  to  build  up  exponentially  as  a function  of  time. 

2)  Following  a thrust,  in-track  position,  and  consequently  the 
ground  track,  rapidly  diverge  from  the  no-thrust  predictions. 

3)  SCF  and  ADC  elements  result  in  comparable  orbit  predictions. 
No  systematic  bias  is  evident. 

4)  The  first  SCF  element  after  each  thrust  is  unreliable,  pro- 
bably because  some  pre- thrust  conditions  are  factored  into  the  orbit 
determination. 

5)  The  most  reliable  positional  information  is,  of  course,  ob- 
tained close  to  the  element  epochs.  A LOKANGL  interpolation  run  using 
all  elements  naturally  matches  this  orbit  closely  (compare  in  Table  1 
against  diagonal  entries).  If  one  element  set  only  is  used  in  each 
thrust  Interval,  the  prediction  gradually  diverges  as  expected. 

6)  A separate  ROPP  test  showed  that  for  an  S3-4  type  satellite 
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TASlE  2.  R7»P  RUNS  WITH  « DC  ELEMENTS 


a 20%  error  In  the  estimate  of  drag  produces  a 0.7°  mean  anomaly  or 
10  second  In- track  timing  error  after  2 days.  Due  to  errors  In  deter- 
mining the  semi-major  axis,  greater  variations  are  likely. 


1.3  Satellite  Thrust  F.ffects  Study 

The  effect  of  different  amounts  and  directions  of  thrusting  on  various 
satellite  orbits  is  an  extensive  subject,  and  only  the  special  case  of 
the  S3-4  satellite  was  explored.  The  nearly  circular  (e».006), 

96.5°  inclination,  low  altitude  orbit  is  maintained  with  argument  of 
perigee  around  120°  by  periodic  thrusts  that  cause  a velocity  Increase 
of  the  order  of  7.5  m/sec,  approximately  in  the  direction  of  the  velo- 
city vector. 

The  effect  of  a 0.1%  increase  in  the  velocity  vector  on  the  argument 
of  perigee  and  the  eccentricity  was  evaluated  at  30°  Intervals  of  mean 
anomaly.  The  results  are  presented  In  Table  3.  In  all  cases  the  mean 
motion  decreased  from  16.28  to  about  16.231  revs/day  and  the  mean  anom- 
aly change  complemented  the  change  in  the  argument  of  perigee. 

This  table  can  be  used  to  confirm  the  effects  of  known  thrusts  at  var- 
ious latitudes,  ascending  and  descending,  and  eventually  to  Incorporate 
simple  thrust  models  Into  LOKANGL.  For  example,  a regular  burn  occurs 
every  three  days  near  0°  latitude  descending,  with  a resulting  0.1% 
Increase  In  velocity.  The  Increase  In  the  argument  of  perigee  and  the 
corresponding  lower  perigee  latitude  makes  up  for  the  secular  *4°/day 
drift  In  the  perigee.  As  another  example,  on  21  March  1978  two  burns 
gave  a 0.07%  Increase  In  velocity  near  60°  latitude  descending  (M«0°), 
followed  by  a 0.01%  velocity  Increase  near  60*  latitude  ascending 
(M«300°).  The  result  was  a minor  change  in  the  argument  of  perigee. 
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Table  3.  Effect  of  0.1%  Velocity  Increase  on  a Satellite  with 
1 * 96.5°,  e * .006,  u * 120*,  and  n * 16.28 


Original 

Mean  Anomaly  # 

New  - Old 

Argument  of  Perigee 

New 

Eccentricity 

0 

+ 0.1 

.0080 

30 

+ 6.7 

.0079 

60 

+ 12.5 

.0073 

90 

+ 16.3 

.0065 

120 

+ 15.9 

.0055 

150 

+ 12.2 

.0046 

180 

- 0.4 

.0040 

210 

- 10.8 

.0043 

240 

- 16.3 

.0051 

270 

- 16.3 

.0062 

300 

- 12.4 

.0071 

330 

- 7.3 

.0077 
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1.4  ROPP  Features 

The  input  facilities  described  in  Reference  3 have  been  modified  over 
the  past  several  years  as  additional  features  have  been  incorporated 
in  program  ROPP.  Drag  optimization  for  two  element* sets,  and  input 
and  output  of  position-velocity  vectors  was  described  in  Reference  2. 
Inclusion  of  geopotential  resonance  terms  was  described  in  Reference  4. 
Figure  6 reflects  the  up-to-date  status  of  input  facilities  and  options 
for  ROPP. 


0**PLATCCOCO,  0,9)  ROPP 

PR3GPAH  POPP (INPUT, OUTPUT, T A PP5 5 IN»UT, TAPE3sOU^°t**TaPE7f tape 2t  ROPr 

1 T A®?  S * OUT PUT  t PUNC*)  ROPF 

RGPF 


f |f<  i£  k «1  O »<  N ^ ^ |0  K H/T'O^MM^U'>flk«ff'O-<N^#l^.0k'0 


u.  ti  uuvkLi.a.OkA.u.uuu.u.u.u.uau.u.u.ua.u.u.u.u.uu.uuuUuuuu.n.u.uuu.u.uUii.  a u u u.  o u.  u. 
&Q.OLaa.u.&&.a.ixa.iiAii(LiL(L&'(i(LiiiXft.aii.o.(LO.iL(i&a.aQ.a.(Laaaa(iQ.aIaaaaQ.auflLii.ii 
oooooooooooooooooooqogooooooooooooooooooooooooooooooo 
ororaro'a'a'nrflrarn'nrara'aroraro'nfnrnrorarft'oro'oraro'arrt'o'o'aro'o'ortt'nf'n'arnrorartfflrrva'aro-'n'a'nrnr 


7 «r  < 

;>  » 4 r 

»t  n o « c 

J HIU 

a j ► 

3 0 0 » 

u.  w a o 


j «/>  < 

► 3 • C 

-4  I Jtfl 

to  > *f  « 

• OU1 

U « (Lit 

Z to.  * • \f\  r 

o OT  I l 

*-*  uiP  ar  O’ 

► X O O ( 

a 3trt  c 

M • Z Z • c « 

(X  U4  I Z U < 

u s Z 3 < « . 

m « a j u o • 

«4i  z jo  r i i 

C • 4U  ••Ml 


X 

- o 

to.  « O 1-4 
Ui  -J  Ui  to 

a.  c r 
x o ~ 

; H Z 
2 to  o o 
w Z H 

KCCh 

to  O -i  — 
to  u »- 
z o U'  » r 


HMO'k 
I I I I 


J J J J 
uuoo 
ogou 


l/l  m 4T. 

ui  * x uj 

to  OOH 

« j«*  or 

o c _l  to  ►- 

MUiO  z 
o ar  to.  ui  ui 

Z to»  to 

poaro 

to  * 

o m a ar 

•u  to 
- -i  m a» 

• to  « r 
to  3 3 in  D 
VI  to  CJ  ^ Z 
JZUZ 

E H T to 

7 «•  z o r 

• I «*  o •* 


: z < o nr 

- «*  o a ui 
r to  ♦ o > 

> J IT*  Z too 

> D • U to  VO 

->  o tr  a to  i 

* to  «V  4« 

. u to  «4  to  to  m ui 

I •*  x <r  o in  to  x • 

3 to  »-4  to  to  CJ  to  ; 

•IT  U ^ or  « MU  l 

3 to  ivi  uj  a or 
’ wo  » novai 
» u o u r to  Ui  o 

t u.  to  to  =>  VO  > to  I 

- o m c z to  o %r  * 


Z to  -i 

3 «*  ■« 

O > > 

x or  or 

• Ui  to 

• to  to  ►- 

to  < Z Z 

Z3  O to  to  l 

a •.  • i 

Z O o to  to  i 

to  r h z r ; 


to  to  Z 
Z < to  X 
to  ti  r •-< 
►-  to  to  r 

or  or  or  y 
a a a z 


OUOOOOOOOUOUOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOUOtJOOOOO 


THIS  PAGA  IS  BEST  QUALITY  PRACTICAL 
(X)|PZ  JJRMISHSI)  XQDD.Q 


37 


nerMOB.fl  ROPP 

SEMIAX*  r FENI-iAJCR  AXIS  (KM,  Or  F.R.  IF  ROPF 

ANULT  *1,1  ROPF 

ECCEN*  F ECCENTRICITY  ROPF 

INCL  • F INCLINATION  (0£GRrE*)  ROPE 


76/76  OPT*t  rT*  6. **65?  07/10/7*  16.0a. 07  P»GF 


Figure  6c.  Program  ROPP  Input  Format 
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1.5  Eclipse  Calculations:  Program  ECLLUH 

\ I ' 

Synchronous  orbit  satellites  have  exhibited  a tendency  to  acquire,  at 

times,  high  levels  of  electric  charge.  The  associated  elevated  voltage 

i,  i 

I i levels  can  cause  breakdown  with  resulting  adverse  effects.  Insulating  * 

materials,  for  example,  can  undergo  Irreversible  damage.  In  addition, 
the  Impulsive  flow  of  current  accompanying  breakdown  can  Induce  exces- 
sive electrical  transients  In  on  board  electronics,  resulting  In  tempo- 
rary or  permanent  malfunctions. 

V,*.  : Y ••  ■ I 

Charge  build-up  occurs  as  a result  of  Interaction  of  the  satellite  with 
Its  plasma  and  solar  radiation  environments.  The  solar  radiation  level 
undergoes  extreme  variations  as  the  satellite  proceeds  Into  and  emerges 
from  conditions  of  eclipsing  by  the  earth.  To  analyze  the  role  of 
varying  Illumination  Intensity  In  the  charging  process,  one  needs  to  be 
able  to  calculate  Illumination  on  the  satellite  as  a function  of  time. 

This,  In  turn,  requires  precise  evaluation,  as  a function  of  time,  of 

' 

the  relative  geometry  of  the  sun,  the  earth,  and  the  satellite.  An 

■ 

accuracy  of  0.1  sec  of  time  Is  the  objective  In  these  calculations. 

Two  approaches  to  these  calculations  are  available  In  terms  of  existing 
software.  The  first  Is  the  routine  ECLIPS  which  Is  used  In  conjunction 
with  ROPP  and  Is  documented  In  Reference  4.  The  second  Is  a routine, 

| 

ZWEI/INTEG,  prepared  by  the  researcher  and  documented  In  Reference  5. 

Neither  Is  Inherently  adapted  to  the  requirement  for  0.1  sec  precision. 

In  addition,  these  programs  perform  distinctly  different  functions. 

) 
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ECLIPS  solves  a transcendental  equation  to  obtain  two  numbers:  the 
times  at  which  the  satellite  enters  and  exits  the  earth's  shadow. 

ZWEI/INTEG  computes,  at  successive  Instants  of  time,  the  percentage  i 

i 

obscuration  of  the  sun  as  It  sets  and  rises  behind  the  earth  limb. 

Since  It  Is  the  UV  component  of  the  solar  radiation  that  Impacts  the  i 

i 

charging/discharging  process,  ZWEI/INTEG  takes  account  of  the  substan- 
tlal  UV  absorption  In  the  earth's  atmosphere.  ECLIPS,  on  the  other 
hand.  Ignores  the  atmosphere  and  calculates  on  the  basis  of  shielding 
by  the  earth's  surface.  The  objective  of  the  present  task  Is  to  Iden- 
tlfy  sources  of  Inaccuracy  and  develop  Improved  calculation  for  eclipse 
phenomena. 

1.5.1  Factors  Affecting  Accuracy 

ZWEI/INTEG  performs  essentially  two  functions.  First  are  the  solar  and 
satellite  ephemerls  calculations  that  establish  the  geometry.  Second 
are  the  absorption/shadowing  calculations  Involving  an  Integration  over 
the  unobstructed  portion  of  the  solar  disk  using  an  atmospheric  absorp- 
tlon  model.  Included  In  this  calculation  Is  the  effect  of  variation  of 
atmospheric  absorption  with  altitude.  The  objective  of  the  work  re- 
ported here  Is  to  upgrade  the  first  calculation  while  leaving  the  sec- 
ond Intact.  Program  ECLLUM  performs  this  function. 

Factors  which  affect  the  accuracy  of  eclipse  calculations  Include: 

1.  Precision  of  calculations  of  satellite  coordinates; 

2.  Oblateness  of  the  earth; 

3.  Solar  ephemerls  accuracy; 

4.  Time  for  which  solar  coordinates  are  evaluated; 
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5.  Assumption  that  the  satellite-sun  line  Is  parallel 
to  the  satellite-earth  line. 

Satellite  orbit  calculations  must  be  performed  to  provide  the  changing 
, location  of  the  satellite  as  It  flies  Into  and  out  of  the  shadow  re- 

gion. ECLLUM  employs  a simple  Keplerlan  orbit  over  the  several  minute 
time  spans  of  Interest  during  which  the  sun  rises  and  sets,  as  observed 
In  the  UV  regime,  at  the  satellite.  Orbital  elements  for  these  cal- 
culatlons  are  obtained  by  using  ROPP  to  evaluate  the  elements  applic- 
able at  eclipse-in  and  -out  times  provided  by  the  subroutine  ECLIPS. 

To  Insure  accuracy,  mean  elements  are  not  used;  rather,  the  osculating 
i (or  Instantaneous)  elements  provided  by  ROPP  and  evaluated  at  the  esti- 

mated eclipse  times  are  employed.  A discussion  of  some  of  the  sources 
of  error  and  the  error  magnitudes  Inherent  In  using  ROPP  predictions 
has  been  presented  above  In  Section  1.2. 

Ignoring  oblateness  of  the  earth  and  treating  It  as  a sphere  of  radius 
equal  to  the  earth's  actual  equatorial  radius  can  Introduce  up  to 
approximately  20  km  error  In  the  location  In  space  of  the  region  of  the 
earth's  surface  responsible  for  eclipsing.  In  a worst  case  situation 
this  could  lead  to  a six  or  seven  second  time  error.  To  avoid  this 
type  of  potential  error,  ECLLUM  takes  specific  account  of  the  earth's 
oblateness. 

I I 

A major  factor  affecting  the  precision  of  the  solar  coordinates  Is  the 

it 

sophistication  of  the  algorithm  employed.  As  discussed  In  Section  2.3 

i ' 
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three  routines  have  been  evaluated.  To  meet  the  requirement  for 
timing  precision,  the  routine  SOLEPH  has  been  selected  for  use  in 
ECLLUM.  It  is  expected  to  contribute  errors  no  greater  than  0.03  sec- 
onds of  time. 

Obviously,  even  the  best  algorithm  will  provide  inaccurate  solar  coord- 
inates If  it  is  not  evaluated  at  the  proper  Instant  of  time.  To  Insure 
accuracy  of  the  solar  coordinates,  ROPP/ECLIPS  is  run  to  yield  approxi- 
mate eclipse-in  and  eclipse-out  times.  Solar  coordinates  are  then 
evaluated  at  these  two  times.  This  Is  actually  the  second  step  of  an 
Iterative  approach,  since  ECLIPS  Itself  must  evaluate  solar  coordinates; 
without  a priori  knowledge  of  eclipse  times,  however,  ECLIPS  calculates 
these  coordinates  using  only  very  approximate  times  Inputted  to  a very 
approximate  algorithm. 

To  simplify  calculation,  the  assumption  Is  sometimes  made  that  all  rays 
from  the  sun  to  the  satellite  and  to  points  on  the  earth  originate  from 
a point-source  at  Infinity  (l.e.,  all  rays  are  parallel).  However,  this 
can  Introduce  a timing  error  of  up  to  0.5  sec.  In  a calculation  Inten- 
ded to  provide  0.1  sec  precision,  this  luxury  cannot  be  afforded,  and  a 
more  exact  calculation  Is  required.  The  approach  employed  In  ECLLUM, 
Including  evaluation  of  the  effect  of  oblateness.  Is  presented  In  the 
following  section. 

1.5.2  Analysis 

The  routine  INTEG,  which  embodies  the  physical  model  used  to  evaluate 

j 

solar  obscuration,  requires  as  Input  the  distance  of  closest  approach, 
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Xm,  between  the  surface  of  the  earth  and  the  line  between  the  satellite 
and  the  center  of  the  sun.  In  addition,  the  angle  between  this  same 
line  and  the  earth-satellite  line.  Is  also  required.  Figure  7 illus- 
trates the  geometry,  i can  be  found  from 

i 

| 

r • R r • (R®  - r) 

cos  i - = , 

|r|m  l*--r||7| 

where  and  r can  be  obtained  from  solar  and  satellite  ephemerldes, 
respectively.  It  Is  convenient  to  express  7 In  the  P,  Q,  W coordinate 
system.  (See  Ref.  4,  p.  24,  and  Ref.  6.) 

| 

r = x P + y 0 

O)  0) 

where 

x * r cos  v 

0) 

y = r sin  v 

a) 

P 

r • 

1 + e cos  v 

p ■ a (l-e)(l+e) 

and  e,  a,  and  v are  the  eccentricity,  sem1*major  axis,  and  true  anomaly 
of  the  satellite's  orbit,  respectively. 
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In  the  earth*-centered  Inertial  system. 


i 


Rx 

• «•,  - 

X 

0) 

Px 

3 
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< 

Qx 

1 

1 

R„ 
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X 

P„ 
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Qu 
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i 

y 

y 

(A) 
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0) 
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>» 
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Qz 
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The  orbital  elements  fl,  u,  and  1 used  to  calculate  F and  TJ,  together 
with  a,  e,  and  v used  to  calculate  r,  are  all  obtained  from  ROPP  as  the 
osculating  elements  at  the  approximate  time  of  eclipsing  determined  by 
the  ROPP  subroutine  ECLIPS.  With  the  value  of  ^ determined,  can  be 
found  from 


X$  * r | si n i\ 

The  earth's  radius,  R£,  at  the  shadowing  point,  Q,  Is  provided  as  Input 
to  the  routine  INTEG,  which  then  calculates 

" XS‘RE 


XM  Is  the  distance  of  closest  approach  between  the  earth's  surface 
(assumed  spherical)  and  the  sate!11te*-sun  line.  To  correct  for  the 
earth's  oblateness,  the  value  of  the  earth's  radius  Inputted  to  INTEG 
should  be  the  value  appropriate  for  the  latitude  of  that  portion  of 
the  earth's  limb  responsible  for  eclipsing  at  the  satellite's  point  of 
entrance  to  or  exit  from  the  shadow  zone;  l.e..  It  Is  not  sufficient 
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simply  to  use  the  earth's  equatorial  radius 


Fig.  8a  Illustrates  the  geometry  of  the  oblate  earth.  It  is  assumed 
that  the  earth  can  be  represented  as  an  oblate  spheroid  with  its 
semi-major  axis  given  by  the  earth's  equatorial  radius,  A£.  All 
meridian  cross  sections  are  ellipses  of  semi-major  axis,  A£  and 
eccentricity  eE.  All  cross  sections  in  planes  parallel  to  the  equa- 
tor are  circles. 


Point  Q in  Fig.  8a  is  assumed  to  be  the  shadowing  point  on  the  earth's 
surface.  The  eclipsing  calculations  require  the  earth's  radius  at 
that  point.  The  approach  to  be  used  is: 

a)  Evaluate  az,  the  z-component  of  point  Q.  (i.e.,  its 
perpendicular  distance  from  the  equatorial  plane.) 

b)  az  defines  the  geocentric  latitude,  d>,  of  Q. 

c)  Knowledge  of  <p  then  permits  calculation  of  the  radial 
distance  from  the  earth's  center  to  Q. 

To  evaluate  az,  observe  from  Fig.  8b  that 
7 = a + 1 
where  111=  Rr 


and  1 Is  that  generator  of  the  cylinder  bounding  the  nominal  (visible) 
shadow  region  which  passes  through  the  satellite's  point  of  entrance 
Into  the  shadow  region. 

Therefore 


a = r - d 
z z z 


Assuming  a point  source  sun  at  Infinite  distance  (an  assumption  that 
is  valid  for  earth  flattening  calculations,  but  not  elsewhere  in  this 
analysis)  it  is  clear  that  the  vector  d is  in  the  reverse  solar  dir- 
ection. Therefore 


az  = rz  + ld 


cos  ~ 


Therefore  az  = rz  + r cos  ^ 


But  from  Keplerian  orbit  analysis. 


1 + e cos  v 


and  therefore. 


az=  P 


(Pz  cos  v + Qz  sin  v + jjj^  cos  ^ 


1 + e cos  v 


From  the  properties  of  ellipses. 


(1  - eE  sin. 


f 

I 

1 

I 

2 

When  this  is  solved  for  cos  <P, 

! 

AE2  t1  - eE2]  “ az2 

2 

COS  <p  = j 

AE2  [1  - eE2]  - az2  eE2 

» 

p 

To  obtain  R^,  the  radius  at  point  Q,  this  value  of  cos  Is  substi- 
tuted Into 


At  that  point,  all  required  Inputs  for  INTEG  are  available: 

X$,  'P,  and  the  value  of  RF  at  point  Q. 

A grid  of  equl -spaced  time  Instants  Is  established  which  spans  the  en- 
tire eclipsing  process  (l.e.,  extending  from  zero  to  total  obscuration). 
At  each  Instant,  the  satellite's  location  Is  evaluated,  using  simple 
Keplerlan  orbital  dynamics  based  on  osculating  elements  evaluated  by 
ROPP  for  the  nominal  times  of  shadow  region  entrance  and  exit.  As 
values  for  X$,  *P , and  R£  are  calculated,  they  are  provided  as  input 
to  subroutine  INTEG,  which  evaluates  the  corresponding  illumination. 

The  result  Is  a plot,  versus  time,  of  variation  of  solar  UV  radiation 
Intensity  Incident  upon  the  satellite  during  periods  of  penetration 
Into  and  emergence  from  the  earth's  shadow. 
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2.  Ionospheric  Research  Support 
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2.  Ionospheric  Research  Support 


Programs  to  describe  environmental  conditions  tnat  are  significant  to 
ionospheric  phenomena  have  been  developed  at  AFGL  and  other  agencies. 
Further  work  to  acquire  and  extend  capabilities  for  describing  proton- 
electron,  geomagnetic,  and  astronomical  parameters  Is  reported  in  this 
section. 
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2.1  Proton-Electron  User  System 

i 

i 

Initiator:  H.  B.  Garrett 

Project  No:  7661  Problem  No:  4913 

* 

r 

User  system  PROTEL  consists  of  a series  of  data  models  describing  the 
Inner  and  outer  radiation  zone  electron  and  proton  environments  and  a 
collection  of  three  interlocking  main  subroutines  and  various  support- 
ing subroutines  which  perform  the  following  functions: 

1.  Subroutine  ORB  generates  orbit  tapes  and  calculates  B-L  values 
suitable  for  Input  to  subroutine  ORP  using  the  Brouwer  orbit  generator 
suited  to  orbits  with  eccentricities  >0.1. 

2.  Subroutine  ORP  calculates  the  average  geomagnetically  trapped  rad- 
iation accumulated  by  an  orbiting  vehicle. 

3.  Subroutine  MODEL  allows  the  user  to  access  any  of  the  current 
trapped  radiation  models  available  from  the  National  Space  Science  Data 

I 

. 

Center. 

4.  Subroutine  PREMOD  selects  the  appropriate  proton  or  Inner  and  outer 
zone  electron  models  from  the  packed  binary  CDC  EDITLIB  file  PREMOD  for 
use  In  subroutines  ORP  and  MODEL. 

5.  Program  PACKBIN  prepares  the  eight  data  models  for  cataloguing  In 
the  EDITLIB  file  PREMOD  from  the  unpacked  BCD  display  format  used  for 
storage  in  the  CDC  UPDATE  program  library  file  PROTEL.  Data  models  are 
packed  and  written  In  binary  form  on  TAPE8. 

Models  now  available  for  use  at  AFGL  Include  inner  zone  electron  model 


i 

I 
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AE5  (AE5MIN),  solar  minimum,  at  epoch  October,  1967;  inner  zone  elec- 
tron model  AE6  (AE6MAX),  solar  maximum,  epoch  1980;  outer  zone  electron 
models  AE4,  solar  minimum  (AE4MIN)  at  epoch  1964;  and  solar  maximum 
(AE4MAX)  at  epoch  1967;  interim  outer  zone  models  AEI-7,  at  epoch  1980 
s . (AEI7HI  and  AEI7L0);  and  2 AP8  proton  models  at  solar  minimum  epoch 

1964  (AP8MIN) , and  solar  maximum  epoch  1970  (AP8MAX).  The  two  AEI-7 

I I 

models  differ  in  that  AEI7HI  favors  Vampola's  fits  to  the  OUI-19  data 
and  AEI7L0  is  more  representative  of  all  of  the  data  sets  available  to 
NSSDC  as  of  November,  1972.* 

2.1.1  Functional  Description 

The  complete  set  of  all  PROTEL  subroutines  and  unpacked  data  models  are 
accessible  in  CDC  UPDATE  format.  CDC  EDITLIB  subroutine  library  PRELIB 
contains  the  ORB,  ORP,  MODEL,  PREMOD  and  PACKBIN  subroutines.  The 
NSSDC  programs  ORB,  ORP  and  MODEL  were  converted  for  use  on  the  CDC 
6600  NOS/BE  system  through  the  conversion  of  IBM  360  double  precision 
variables  to  single  precision  60-bit  CDC  word  storage  convention  and 
required  character  conversions.  Duplication  of  subroutines  used  by 
more  than  one  of  the  original  programs  was  eliminated  by  the  EDITLIB 
file  construction.  Unnecessary  Input  parameters  were  removed  and  other 
program  inefficiencies  were  corrected  in  order  to  minimize  core  re- 
quirements and  program  execution  time.  Data  models  were  packed  In 
t binary  format  with  three  data  words  per  60  bit  CDC  word,  and  statement 

function  MAP  was  designed  to  unpack  the  appropriate  word  and  select  the 
relevant  bits  as  needed  during  execution.  In  this  way,  run-time  model 
j • storage  space  was  reduced  by  one- third.  NSSDC  models  were  provided  In 

‘ ■ one  of  two  forms:  as  Fortran  IV  BLOCKDATA  packages  or  as  BCD  formatted 

flies  of  nine  words  per  record.  Consistency  was  obtained  through  the 
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elimination  of  BLOCKDATA  format  and  the  use  of  a call  to  subroutine 
PREMOD  to  load  the  requested  proton  model  or  electron  models  Into  the 
designated  storage  area.  Either  one  proton  model  or  an  Inner  and  outer 
electron  model  pair  may  be  called  by  PREMOD. 

2.1.2  Usage  of  the  PROTEL  System 

Figure  1 (a  A b)  summarizes  the  rules  for  usage  of  the  Proton*Electron 
User  System.  Subroutine  call  parameters  consist  of  two  types:  those 
which  declare  input  to  the  subroutine  (subject  to  definition  by  the 
user);  and  those  which  control  the  nature  and  extent  of  the  subrout* 
ine's  output. 

There  exist  certain  Increment  and  range  restrictions  for  the  Inner 
energy  (E)  and  logarithmic  equatorial  magnetic  flux  (L)  values  allowed 
by  subroutines  ORP  and  MODEL.  The  Incremental  restrictions  are  as 
follows: 

Protons: 

for  E < 1 MeV,  AE  2 250  keV 
for  1 S E/(MeV)  £ 20,  AE  2.  MeV 
for  20  < E/(MeV)  i 50,  AE  2 5 MeV 
for  E > 50  MeV,  AE  2.  10  MeV 
Inner  Zone  Electrons: 

for  E < 100  keV,  AE  2 50  keV 
for  100  £ E/(keV)  £ 250,  AE  > 100  keV 
for  E > 250  keV,  AE  > 200  keV 
Outer  Zone  Electrons: 

same  as  inner  zone  electrons,  except 
E > 4 MeV,  AE  2. 100  keV 
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proton-electron  user  system 


USER  CO'JTKO-  CARDSl 
_ ATTACH, PRlLID,  ID  = LOGICON. 

ATTACH,TAPEc,  ,PPtMOD,IO  = LDGICOH. 

_ i iboary,?relir. 


NOTE  I (i)  1APE13  CONTAINS  ORBIT  ANU  B &L  PARAMETERS,  AND  MAY  3E 
CATALOGUED  IF  ORP  IS  RUN  LATtR.  THIS  BINARY  FILE  IS 
GENERATED  BY  ORB  FO=>  USE  BY  ORP.  NOT  REQUIPED  BY  MJO^L. 


(?)  T APE  fc  CONTAINS  ALL  PROTON  ANU  ELECTRON  MODELS  INCLJD1NG 
APS  M AX , AEI7HI, AF  ITLO, AE5MIN ,aE6MAX»  t APJ  MIN • THISFILF  IS 
REQUIRED  ONLY  PY  ORP  OR  MODEL. 

(3)  PRINTED  OUTPUT  I?  TAPF6. 


FIELD  LENGTH?  RFQUIREJ  TO  EXFCUTEl 
OrB  - 46K8 

__ORP  - 46< i (INCLUOFS  PREMOO) 

MODEL  - *6Ko  (INCLUDES  PREMOO) 
ORB, PREMOO, ORP, MOUcL  RUN  TOGETHER 


1 i6K8 


CALL 


IBMVV* 

Or\3(SLC;>PN,OSEC,  XN,  XE,XI  ,X  M , X P , X R,  I PR , LI  M BL ) 
StoSPN  ORBIT  GENERATION  SPAN  IN  SECS. 


XHiS  PAGE  IS  BEST  QUALITY TRA.WIC1&1 
2BQji|  O0|PY  tfURNISHED  TO  DDC  ' 


JiEC 


XN 
_XE 
XI 
_XM 
XP 
_XR 
I PR 

_LI_MB_. 


ORBIT  GENERATION  SPAN  IN  SECS. 

(NSSOC  RECOMMENDS  20+  ORBITS  FOri.GOOQ  B-L  SPACE  COVERAGE) 
TIME  INTERVAL  BETWEEN  POINTS  OF  GENERATED"  ORBIT 
(IF  OStC=  0 • , OF FAJL T DSEC=2 0 . ♦ . 0 331* ALT ( KM ) > 

MEAN  MOTION(RcVS/DAV)  Ok  PEkIGEE  ALTITUOE(KM) 

ECCENTRICITY  OJ.O) 

INCLINATION  (DEG)  Ou.O) 

MEAN  ANOMALY  ( OEG) 

ARGUMENT  OF  PEPIGEE  (OEG) 

RIGHT  ASCENSION  ( DFG ) 


SUPPRESSES 
LIMITS  B-L- 


PRINT-OUT}  =1  PRINT-OUT 
VALUES  ON  TAPE10  ? =1  NO 


REQUESTED 

LIMITATION  ON.  VALUES 


Figure  la.  Proton-Electron  User  System 
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CALL  RR“HOD( 6 HI NR MOD » 6H0U T 100) 

WHERE.  INRMOJ  OPTIONS  ARE 

OJIMOO  OPTIONS  ARP 

ONEHOJ  . uPTIONS  _4R£L 

. 0 AD S SPE CIFIEO  MODELS 


OR  CALL  PREMD D (6 HONEMOD* 1H  ) 
. AE5.M  IN* AE6MAX 
AET7HI > AEI7L0 

APttMIN,AP8MAX  . ... 

FROM  TAFE8. 


CALL  ORP  (TABLOG, Nt,tT, ISFIP,E, IORB) 

AQijJ  IA.L-T  A BLQG  ( 5 ) - rACH  .TRUE.  OR  .FALSE. 

1 INTERMEDIATE  PRINTOUT  TABLE 

2 L-BAND  SJMMAPY  TABLE 

3 INTEGRATED  FLUX  TABLE 

j».  INTENSITY  SUMMARY  TABLE 

5 °E AK  CLUX  TABLE 

OiMENSION  E (NE) , I0R9 (2) 

NE  NUMBER  OF_EN"RGv  VALUES  '£*  .{<3L)__ 

ET  THRESHOLD  ENERGY  FOP  TABLES  1 & 5 

ISKIP  READ  1 POINT  LSI  ISKIP  FROM  JABEM-. 

E ENERGY  THP^HDLD  ARP  AY- ME  V (ASCENDING 

. I OF  3 INDEX  NUMBERS  FOR  FI  PST  4 l Ail  ORBITS 


CALL  MODEL  ( NE , NL ,1 OIFF , ID"F , E , XL , Bu 1 ,B0 2 , NOEL 8) 
OIMENSIOT  t ( NE  ) * XL ( NL ) , 3jT( N L)  , BO 2 (NL)  , NDeTbTnL  ) 


ORDER) 


NE  NUMBER  OP  ENERGf  VALUES  »E*  (<10) 

IL  _.  .NUMBER  OF_.l^_VAj^lES  ..(.<ljU) ■ 

I Dir F 0=  INTEGRAL  FLUX  UUTPUTJ1  = AVG.  DIPFERENTIAl  rLUX  OUTPUT 

_ 2=  BU  TM  PUTPUIS _ 

E ENERGY  TMRPSHOLO  ARRAY  - MEV  (ASCENDING  ORDER) 

_ .XL L VALUE  ARRAY 

I 0Er  DEFAULT  ’5  TD  30  LINEAR  B/90  INCREMENTS  TD  CUTOFF 

.1=1  NPUT^i1,3.c.SxmLB_,)fAULES- 

NDTEf  FOR  IDEF=1,  ONF  SET  BELOW  PER  L VALUE- 

BQ1  LOWER.  ^IflITLJIF_9/_aO. 

B 0 2 UPPER  LIMIT  OF  B/BO 

M 


f 

I 

» 

Restrictions  on  the  range  of  values  Include  the  following: 
for  AE5MIN,  L i 1.2  R 

e l 

for  Inner  zone  electron  models,  E i 4.5  MeV 

Although  usage  of  inner  and  outer  electron  models  in  simultaneous  runs 

i 

are  the  most  physically  accurate,  they  are  not  required  by  the  PROTEL  > 

software  system.  In  any  case,  simultaneous  usage  of  inner  and  outer 
zone  electron  models  requires  that  inner  zone  electron  parameter  res- 
trictions must  be  followed.  Failure  to  restrict  parameters  to  the 
above  specifications  results  in  abnormal  program  termination.  In  sub- 
routine ORB,  two  range  restrictions  exists.  Variables  XI  and  XE  may 
not  equal  zero,  and  In  the  case  of  zero-inclination  or  zero-eccentri- 
city orbits,  values  In  the  range  of  10*"9  < XI,  XE  < 10^7  should  be 
entered. 

Subroutine  ORB  allows  the  user  the  choice  of  a point-by-point  tabular 
printout  of  the  generated  orbit  parameters  (Fig.  2)  and  an  abbreviated 
printout  in  which  the  first  and  last  few  points  only  are  displayed. 

Subroutine  ORP  has  the  capability  to  generate  the  foil wing  tabular  out- 
put: 

Intermediate  Printout  (Fig.  3)  - a point  by  point  table  of  the 
omnidirectional  Integral  flux  at  each  point  of  the  orbit  for  a 
given  threshold  energy. 

L-Band  Summary  (Fig.  4)  - a summary  of  the  omnidirectional  part- 
icle flux  (partlcles/cm^-day)  accumulated  In  arbitrary  energy  and 
L bands. 
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Figure  2.  Printout  for  ORB 
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Integrated  Flux  (Fig.  5)  - a summary  of  the  Integrated  flux 
accumulated  In  arbitrary  energy  bands. 

Intensity  Summary  (Fig.  6)  - a summary  of  the  omnidirectional  i 

I 

particle  flux  accumulated  in  arbitrary  energy  and  intensity  bands. 

Peak  Flux  per  Orbit  (Fig.  7)  - a table  of  peak  omnidirectional  i 

f 

integral  flux  encountered  for  each  revolution  for  a given  energy 
threshold. 


Subroutine  MODEL  allows  the  user  the  choice  of  integral  (Fig.  8)  and 
differential  (Fig.  9)  flux  tables  as  output.  Subroutine  PREMOD  prints 
an  informative  message  which  confirms  the  model (s)  loaded  and 
describes  any  input  errors  which  m^  have  occurred. 


Creation  of  packed  binary  file  PREMOD,  which  holds  the  trapped  radia- 
tion models.  Involves  the  use  of  program  PACKBIN  and  the  UPDATE  program 
library  PROTEL.  The  following  NOS/BE  control  cards  are  required  to 
execute  this  mode  of  the  system: 

ATTACH,  OLDPL,  PROTELPL,  ID=L0SCAG. 

REQUEST,  TAPE8,  *PF. 

UPDATE,  C=TAPE7. 

UPDATE. 

FTN,  SL,  R,  I. 

LGO, PL *77777. 

CATALOG,  TAPE8,  PREMOD,  ID=L0GIC0N. 

7/8/9 

*C  AP8MAX.AE4MIN 

7/8/9 

*C  PACKBIN 

7/8/9 

6/7/B/9 
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Figure  6.  Printout  of  Distribution  of  Particle  Flux 

Within  Energy  Bands  and  Intensity  Bands  (ORP) 
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1. 

INNS?  ZONE  ELECTRON  MODEL  »E6  SO.*?  MAX  PROJECTED 
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Figure  9.  Printout  of  Differential  Flux  (M0DE1) 


Program  PACKBIN  supplies  a printout  of  the  contents  of  the  packed  bin- 
ary file. 
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2.2  Geomagnetic-Geographic  Coordinate  Conversion 
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Initiator:  J.  Buchau 

Project  No:  4643  Problem  No:  4965 

Geographic  to  geomagnetic  coordinate  conversion  and  its  inversion 
require  the  availability  of  a one-to-one  transformation  model  for 
latitude  and  longitude.  The  centered  dipole  model  with  the  prime  geo- 
magnetic meridian  containing  the  south  geographic  pole  is  the  simplest 
approximation  for  this  purpose,  where  the  location  of  the  dipole  poles 
are  prescribed  for  the  epoch  by  any  selected  field  model.  Since  the 
reference  longitude  and  latitude  are  hereby  established,  transformation 
from  geocentric  to  geomagnetic  coordinates  follows  by  application  of 
spherical  trigonometry. 

5 6 

Hakura  , Gustafsson  , and  others  have  introduced  "corrected"  geomag- 
netic coordinates  to  account  for  the  asymmetric  field  lines  represented 
by  the  higher  order  spherical  harmonic  terms  in  any  field  model  such  as 
GSFC  12/66,  IGRF  (1975),  etc.  The  correction  to  the  dipole  coordinates 
is  given  by  the  displacement,  due  to  the  higher  order  terms,  of  a field 
line  starting  In  the  geomagnetic  equatorial  plane  and  terminating  at 
the  altitude  at  which  the  coordinates  are  desired.  Gustafsson  tabul- 
ates results  for  the  earth's  surface  using  the  GSFC  field  12/66  for 
epoch  1965,  based  on  an  elaborate  analytical  derivation.  A purpose  In 
his  approach  was  to  evaluate  the  contribution  from  additional  higher 
order  terms,  but  this  procedure  could  not  be  readily  Implemented  at 
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AFGL  for  use  with  recent  field  models.  A computerized  technique  for 
obtaining  a constant  altitude  geographic-geomagnetic  coordinate  conver- 
sion based  on  IGRF  (1975)  is  however  worth  Investigating. 

Two  alternate  programs,  MAGCON  and  MAGFLN,  were  implemented  to  tabulate 

I 

geographic  contours  of  constant  geomagnetic  latitude  and  longitude  by 
examining  the  field  line  models.  Both  programs  make  use  of  the  MGFLD2 
magnetic  field  package  which  was  obtained  from  the  SUA  program  library, 
and  which  is  an  updated  version  of  the  package  described  in  Reference  7. 
MGFLD2  calculates  the  L shell  value  and  the  magnetic  field  components 
for  a specified  field  model,  location  and  epoch,  using  either  the 
SHELLG/FELDG  or  the  slower  but  equivalent  INVAR/FIELDG  system  of  rout- 
ines. In  addition,  field  line  tracing  and  conjugate  point  calculations 
are  available  through  subroutine  LINTRA. 

Program  MAGCON  tabulates  geographic  latitude  vs.  magnetic  dip  latitude 
at  5°  geographic  longitude  Intervals.  Since  the  magnetic  dip  latitude 

I 

(DL)  Is  given  by 

tan  (DL)  = 0.5  (£) 

where  Z Is  the  vertical  and  H Is  the  horizontal  magnetic  field  compon* 
ent,  the  geographic  latitude  Is  Iteratively  obtained.  Figure  10  shows 
sample  printout  using  IGRF  (1975)  updated  to  1978  for  an  altitude  of 
100  km. 


I 
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Program  MAGFLN  tabulates  the  geographic  longitude  trace  of  field  lines 
starting  at  40°N  magnetic  dip  latitude  and  100  km  altitude,  vs.  geo- 
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Geographic  Latitude  vs.  Magnetic 
Altitude  - Sample  MAGCON  Printout 
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detic  latitude  at  5°  geographic  longitude  intervals.  Figure  11  shows  a 
sample  print-out.  Subroutine  LINTRA  was  modified  by  an  entry  point 

i 

CONTRA  which  returns  to  the  calling  program  at  any  specified  geodetic 

r 

latitude  GDLATF,  but  resumes  field  line  tracing  from  this  point  on  a 

) 

, subsequent  call.  If  GDLATF  > 90°,  a normal  conjugate  point  calcula-  t 

1 

tion  terminating  at  a specified  altitude  is  performed. 
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Figure  lib.  Geographic  Longitude  Trace  of  Field  Line  from  Magnetic  Dip  Latitude  40°N,  100  km  Altitude 
Sample  MAGFLN  Printout 


2.3  Precise  Solar  Ephemeris  Routine 
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Initiator:  E.  Robinson 

Project  No:  0001  Problem  No:  4517 

Routines  at  three  levels  of  speed,  core  memory,  and  accuracy  are  avail- 

Q 

able  for  various  solar  and  lunar  ephemeris  requirements.  SOLUN  ap- 
proximates these  ephemerides  using  about  ten  coefficients  for  the  sun. 

A limited  test  showed  that  solar  ephemeris  accuracy  is  75"  arc  in  right 
ascension  and  10"  arc  in  declination  when  compared  to  the  U.S.  Naval 
Observatory  Ephemeris.  This  routine  is  satisfactory  for  many  experi- 
mental data  correlation  requirements. 

Program  SOLLUN  is  a higher  accuracy,  self-contained  package  that  was 
obtained  from  the  SUA  program  library  (Project  0001;  Problems  1131, 
1461),  and  which  is  used  to  provide  most  of  the  print-outs  or  plots 
requested  for  various  observing  stations.  The  solar  ephemeris  is  based 
on  about  thirty  coefficients  and  in  a comparative  test  gave  accuracies 
of  15"  arc  in  right  ascension  and  7"  arc  in  declination. 

8 

Subroutine  SLEPHEM  was  developed  earlier  for  the  highest  accuracy  re- 
quirements at  AFGL,  such  as  eclipse  time  calculations.  For  this  pur- 
pose, celestial  latitude  and  longitude  calculations  are  carried  out 
Including  all  periodic  coefficients  that  are  greater  than  .025"  arc. 

The  sun-moon  coefficients  comprise  789  argument-index  sets  and  are 
stored  on  permanent  file  for  execution  of  SLEPHEM.  A comparative  test 
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showed  right  ascension  and  declination  accurate  to  0.5"  arc.  For 
future  high  accuracy  satellite  eclipsing  calculations,  the  approxi- 
mately sixty  solar  coefficient-index  sets  were  incorporated  into  a 

I 

self-contained  subroutine  SOLEPH.  Right  ascension,  declination,  radius 

vector  and  semi-diameter  are  calculated  for  any  universal  time,  with 

I 

specified  or  default  correction  for  ephemeris  time.  Optional  print-out 


for  checking  against  the  American  Ephemeris  and  Nautical  Almanac  veri- 
fies and  calibrates  the  calculations. 
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3.  Atmospheric  Density  Models 


Initiator:  K.  S.  W.  Champion 


Project  No:  6690 


Problem  No:  4867 


The  continuing  need  for  accuracy  in  satellite  tracking  and  ephemerides 
prediction  in  high  atmospheric  drag  situations  results  in  the  need  for 
improved  modeling  of  the  complex  lower  thermospheric  density  variations. 
A case  in  point  concerns  the  Skylab,  orbiting  space  laboratory,  whose 
orbit  is  decaying  at  a faster  rate  than  predicted  from  density  models 
based  on  the  previous  solar  activity  cycle.  This  clearly  points  up  the 
need  for  better  understanding  of  interactions  of  solar  radiation  with 
the  neutral  atmosphere.  Geomagnetic  activity  remains  the  dominating 
factor  in  short  term  density  variations  and  orbit  determination. 

Recent  in-situ  accelerometer  and  mass  spectrometer  data  provide  fuel 
for  these  studies.  Support  for  this  work  has  taken  three  forms: 

1)  Evaluation  of  recent  density  models  for  ephemeris  prediction 
accuracy. 

2)  Assistance  in  development  of  empirical  density  models  by 
adaptation  of  multiple  linear  regression  software. 

3)  Analysis  of  high  quality  Doppler  beacon  satellite  tracking 
data,  with  program  CELEST,  to  obtain  new  satellite  drag  density  data 
for  new  model  development  and  checking  of  accelerometer  data. 

3.1  Recent  Density  Model  Evaluation 

AFGL  has  an  on-going  program  of  evaluating  atmospheric  density  models 
for  their  satellite  ephemeris  prediction  capabilities.  The  primary 
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software  vehicles  for  this  are  the  CADNIP/BADMEP  package  and  program 
2 

CELEST.  Both  systems  determine  a satellite  trajectory,  using  speci- 
fied models,  by  least-squares  fitting  of  tracking  observations  with 
the  resulting  trajectory  propagated  through  numerical  integration 
into  a prediction  realm.  When  the  tracking  observations  consists  of 
positional  information  (i.e.,  range,  azimuth,  elevation),  as  is  the 
case  with  the  data  handled  by  CADNIP/BADMEP,  then  a direct  comparison 
of  the  predicted  trajectory  with  the  observations  is  possible.  In 
the  case  of  Doppler  data  handled  by  CELEST,  it  is  more  convenient  to 
make  a second  fit  over  the  prediction  realm,  for  comparison  with  the 
predicted  trajectory.  Program  TRDIFF  has  been  developed  for  computing 
and  outputting  differences  between  CELEST -determined  trajectories 
resulting  from  the  use  of  various  models.  As  new  models  are  developed 
they  may  be  added  to,  or  replace,  existing  models  in  the  orbit- determ- 
ination programs.  Recently  this  approach  has  been  applied  to  evalu- 

3 

ation  of  geopotential  models. 

3.1.1  Additions  to  Existing  Programs 

4 

Among  recently  developed  density  models,  the  Jacchia  1977 (J 77),  the 
Mass  Spectrometer  and  Incoherent  Scatter  (MSIS)  model  of  Hedin,  et  al.,! 
and  the  USSR  Apollo-Soyuz  model  have  been  selected  for  Incorporation 
into  the  orbit  determination  programs.  Subroutine  decks  provided  to 
AFGL  by  developing  or  using  agencies  simplified  this  task.  Figure  1 
shows  a schematic  diagram  of  the  Interphasing  developed  In  program 
CELEST.  Subroutine  DENS  had  been  adopted  from  CADNIP  earlier  In 
adding  the  Jacchia  1964  and  other  densities  to  CELEST.7  Hence  the 
Interphasing  for  CADNIP/BADMEP  is  the  same  except  for  the  integration 
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routines  leading  to  the  calls  to  DENS.  Subroutine  READMOD  Is  called 
for  Initialization  of  the  J77  model.  Including  a CDC  Fortran  extended- 
compatible  random  access  mass  storage  file.  The  data  for  this  file  was 
supplied  by  the  developing  agency,  SAO.  Subroutine  APFTMP,  already  In 
existence  for  reading  and  storing  solar  and  geophysical  activity  data, 
was  modified  along  with  DENS  to  provide  the  necessary  Inputs  to  ATMDEN, 
RHO  and  GTS3.  Some  of  the  more  significant  of  these  are: 

Time  (Modified  Julian  Days  and  Fraction)  - Input  to  DENS 

X,  Y,  Z:  Satellite's  ECI  cartesian  coordinates  - input  to  DENS 

DEC:  Satellite  declination  = TAN-1  (Z/>/x2  + Y2) 

ALPHA:  Satellite  right  ascension  = TAN”1  ( Y/X ) 

H:  Satellite  height  = R - RE  [1  - F * SIN2  (DEC)] 
where  R = Jt2  + Y2  + Z2 

RE  = Earth's  equatorial  radius 
F = Earth's  flattening  ratio 
THETAG:  Right  ascension  of  Greenwich 

= THETO  + (TIME-TO)  * THETDT,  with  THETO,  TO,  and 
THETDT  Input  to  DENS  In  initialization  call 
XLONG:  Satellite  East  longitude  = ALPHA  - THETAG 
GLAT:  Satellite  geodetic  latitude 

= TAN-1  | TAN (DEC)  * [1.  + 2 RE  * F/(RE  +H)]} 

GMLAT:  Satellite  geomagnetic  latitude 

• SIN*1  [.9792  S IN ( DEC ) + .2028  COS(DEC)  C0S(XL0NG  - 291°)] 
RSUN , DSUN:  Solar  right  ascension  and  declination,  - obtained  from 
already  existing  subroutine  SOL 
Solar  Obliquity  * 23.44° 
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Sol«r  Distance  * 1 A.U. 
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Solar  and  Geomagnetic  Activity  - From  APFTMP 

3.1.2  Separate  Programs 

Separate  programs,  SOVTEST  and  JAC77,  have  been  produced  as  requested 
for  the  USSR  and  J77  models.  The  Inputs  Include  time  (calendar  date, 
hour,  minute),  height,  west  longitude,  latitude,  solar  flux  (smoothed 
and  instantaneous)  and  geomagnetic  activity  index.  These  are  input 
either  interactively,  using  INTERCOM,  or  non-interactively,  one  card  per 
computation  desired.  In  the  former  case  INTERCOM  types  a header  to  cue 
the  user  for  typing  in  the  data.  From  these  inputs  the  quantities  re- 
quired for  the  subroutines  are  computed  as  described  in  section  3.1.1 
except  that  solar  position  is  obtained  as  described  in  Reference  8. 

3.1.3  Ephemeris  Prediction  Accuracy  Evaluations 

Using  program  CELEST,  updated  as  described  in  section  3.1.1,  satellite 
ephemeris  prediction  accuracy  evaluations  may  be  carried  out  for  the 
various  models.  For  each  density  model  and  time  span,  two  CELEST  orbit 
determination  runs  are  required:  one  for  the  initial  fit  span  termin- 
ated by  integration  of  trajectory  into  the  prediction  period,  and  a 
second  run  determining  a second  trajectory  from  a fit  of  data  over  this 
prediction  period.  The  resulting  trajectories  can  then  be  compared. 

A comparison  of  5 models  thus  requires  10  CELEST  runs  per  case.  To 
streamline  the  procedure  the  system  depicted  in  Figure  2 has  been  devel- 
oped. Program  CELUPF  reads  a set  of  run  parameters  from  cards  (usually 
only  2 cards  are  required).  These  include  parameters  of  existing  pre- 
processed  observation  permanent  file,  if  any,  parameters  of  output 
permanent  file,  day  number  to  be  processed,  model  indicator,  indicators 
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for  any  input  or  output  initial  conditions  permanent  files.  These  are 

combined  with  inputs  from  the  "Daily  Parameters"  file,  which  are  appro-  i 

l 

priate  to  the  day  number  to  be  processed  (daily  polar  motion,  time 

span  to  be  processed,  initial  conditions,  etc.)  to  produce  an  input  J 

I i 

stream  for  the  CDC  UPDATE  facility.  Using  this  with  a standard  CELEST 

I 

input  deck,  in  program  library  format,  as  OLDPL,  UPDATE  creates  the 
BCD  file  COMPILE  which  contains  the  desired  CELEST  run  input  deck 
suitable  for  routing  to  the  input  queue.  The  resulting  run  produces, 
as  its  main  output,  a trajectory  file,  plus  optional  updated  initial 
conditions  and/or  preprocessed  observation  files. 

I When  the  required  trajectory  files  (2  for  each  model)  have  been  pro- 

duced they  may  be  input  to  TRDIFF  to  produce  a printout  of  prediction 

I errors  of  various  components  (radial,  in-track,  cross-track)  for  the 

different  models  for  easy  comparison. 

: ' i 


3.2  Density  Modeling  - Stepwise  Multiple  Regression 

A stepwise  multiple  regression  program  to  produce  empirical  density 

models  from  experimental  data  has  been  developed  from  the  IBM 

. Q 

Scientific  Subroutine  Package  Program  STEPR.  This  program  fits  the 

density  data  to  the  functional  form 

m 

p = a0  + £ aiFi 
i=l 

where  the  a^  are  the  adjustable  coefficients  to  be  determined  and  the 
variables  are  functions  of  various  correlative  parameters  such  as 
solar  flux,  geomagnetic  activity,  etc. 

Models  are  constructed  from  successively  larger  subsets  of  the  vari- 
ables F. , in  each  step  adding  the  variable  which  most  reduces  the 
residuals  between  model  and  data.  After  each  step,  the  latest  model 
is  output  with  statistics  and  the  procedure  terminated  or  continued 
according  to  the  percentage  reduction  in  residuals  in  this  step. 

Figure  3 shows  a block  diagram  of  the  program.  The  major  addition 
is  the  user  supplied  routine  PREPRO  to  generate  the  observation  file 
required  by  DATA.  This  is  called  only  if  the  number  of  observations 
specified  on  the  parameter  input  card  read  by  STEPR  Is  not  positive; 
l.e..  If  the  required  data  file  already  exists,  the  call  to  this 
routine  can  be  skipped.  The  actual  number  of  observations  is  re- 


f 

| 

i 

r 


i 


i 

t 


90 


Parameter 

and 

Selection  Cards 


Figure  3.  Multiple  Regression  Program  STEPR 


turned  by  this  routine  to  STEPR  for  use  by  the  other  structures. 
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The  version  which  currently  exists  is  capable,  with  minimal  modifica- 
tion, of  handling  data  from  satellite  drag  analysis,  accelerometers, 
or  mass  spectrometers.  This  must  exist  on  a file  containing  for  each 
observation:  the  date,  or  modified  Julian  day,  time,  latitude,  longi- 
tude, height  and  density.  Solar  and  geomagnetic  activity  data  are 
assumed  to  exist  on  a file  in  format  accepted  by  CADNIP.*  The  fol- 
lowing variables  are  currently  output 

F10  ^ Solar  flux  average  over  4 solar  rotations 

<KP>12  KP  average  over  previous  12  hours  (Geomagnetic 
activity  index 


semiannual  terms:  sin4  ira,  cos4  ira,  sin4  na  cos2  r a, 
where  a = ( t-36204 ) /365. 2422 

t = time  in  modified  Julian  days 

COS  2 ( X + 30°),  X = West  Longitude 

^10  7 ~ F10  7’  F10  7 = solar  flux  for  previous  day 

<KP>g  KP  average  over  previous  6 hours 

KP(t-3hr) 

Density 

With  minimal  effort  this  set  can  be  modified,  as  is  usually  the  practice. 


92 


« 

I 

1 

? 

3.3  Accelerometer  Data  Support 

Traditional  satellite  orbit  determinations  by  CADNIP  and  CELEST 
provide  useful  checks  or  comparison  for  accelerometer  measurements 
made  on  board  the  same  satellite,  in  addition  to  providing  drag 
density  data.  If  high  quality  Doppler  Beacon  observations  are 
available,  CELEST  is  particularly  useful  because  of  the  capability 
to  subdivide  the  drag  into  4-6  hour  segments. 

Data  has  been  requested  by  AFGL  and  received  from  Defense  Mapping 
Agency  (DMA)  for  satellite  DB-16,  including  the  following: 

1)  Preprocessed  observation  files  in  BCD  format,  for  specific 
time  periods 

2)  Punch  card  deck  containing  data  for  Sun-Moon/coordinate 
transformation  file 

3)  Sample  output  of  CELEST  run 

4)  Thrust  start-stop  times  and  values  (TWX  reports). 

Data  was  transmitted  by  BCD  to  circumvent  Incompatibilities  between 
DMA's  Uni  vac  1100  and  AFGL's  CDC  6600.  This  data  has  been  success- 
fully converted  to  files  needed  by  CELEST  with  programs  DMABIN  (for 
the  preprocessed  observations)  and  SNMNBN  (for  Sun-Moon/transform- 
atlon  data).  Because  preprocessed  data  has  been  supplied  It  will 
not  be  necessary  to  enter  CELEST' s preprocessor  module,  as  was  done 
for  satellites  DB-7  through  DB-9.  Thus  the  rather  large  station 
data  file  Is  not  needed. 

CELEST  is  now  totally  operational  with  this  data. 
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4.0  Introduction 
I 


Program  SSIE  Is  the  main  program  of  a system  to  process  special  Ion 
and  electron  sensor  data  from  a series  of  polar  orbiting  satellites 
of  the  DMSP.  The  program  was  developd  at  AFGL  for  integration  into 
the  real  time  satellite  data  processing  system  at  Air  Weather  Service, 
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Offutt  Air  Force  Base,  Nebraska.  In  the  development  it  was  necessary 
to  simulate  the  AWS  UNIVAC  1110  System  on  the  AFGL  CDC  6600.  The 
structures  of  the  input  and  output  files  for  the  SSIE  program  have 
been  described  In  Reference  1.  This  report  describes  the  Integration 
of  SSIE  Into  the  AWS  System,  and  Into  the  AFGL  operation.  At  AFGL 
the  AWS  system  Is  simulated,  and  In  addition  several  programs  have 
been  developed  to  display  the  data  In  various  ways.  Since  Reference  1 
was  written,  several  changes  have  been  made  to  the  output  file.  These 
changes  are  described  In  Section  4.3. 
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4.1  AWS/AF6L  System 

1 

1.  AWS  System 

1 r 

A schematic  for  the  AWS  System  appears  as  Figure  1.  Currents  generated 

by  electrons  and  positive  ions  are  recorded  from  their  respective  sen*-  t 

sors  on  magnetic  tape.  As  the  satellite  passes  over  ground  stations, 
the  data  is  telemetered  as  the  recorder  rewinds  (hence  backward  in 
time)  to  the  station.  The  ground  station  then  retransmits  the  data 
from  that  readout  (usually  one  orbit)  to  AWS  where  it  becomes  prepro- 
cessor program  input.  The  preprocessor  program  adds  ephemeris  Inform- 
i 

atlon  to  each  minute  of  data  creating  one  minute  data  sets  which  are 
then  placed  In  the  raw  data  circular  file,  which  Is  the  Input  file  for 
the  SSIE  program. 

In  the  pre-operatlonal  phase  the  pre-processed  file  at  Global  Weather 
Central  has  been  dumped  to  magnetic  tape  and  shipped  to  AFGL  to  save 
for  future  studies.  In  the  operational  phase,  SSIE  outputs  the  pro- 
cessed data  to  a circular  file  which  Is  used  for  GWC  modelling  pro- 
grams. Periodically  the  processed  file  Is  dumped  to  tape  and  sent  to 
NOAA  for  archiving.  From  NOAA,  the  data  Is  available  to  a variety  of 
users.  Including  AFGL. 

2.  AFGL  System 

Figure  2 Is  a schematic  for  the  system  at  AFGL.  During  the  pre-oper- 
atlonal phase.  Input  to  the  AFGL  system  has  been  the  tape  copies  of  the 
preprocessed  AWS  circular  file.  When  processed  files  are  available 
from  NOAA,  they  may  be  used  as  Input  after  minor  modifications  to  the 
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Figure  1.  Schematic  of  GWC  SSIE  Data  Flow 
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Figure  2.  Schematic  of  SSIE  Processing  at  AFGL 


100 


f 


) 

I 

mass  storage  file  creation  program.  Figure  3 is  a flow  chart  for  the 
mass  storage  creation  program  GWRIT.  In  GWRIT,  L counts  the  number  of 

I 

1120  UNIVAC  word  records  read,  and  M is  the  count  of  the  112  word  ad- 
dressable blocks  written  on  CDC  mass  storage.  GWRIT  creates  the  input 

to  the  AFGL  version  of  SSIE  on  a private  disk  pack.  The  output  (pro-  , 

cessed)  file  also  resides  on  the  private  disk  pack. 

’ I 

: 

In  the  first  instrument  to  carry  the  ion  and  electron  sensors,  noise 
in  the  monitor  word  amplification  circuit  caused  Irregularities  in  the 
monitor  words  during  the  electron  and  ion  sweeps.  In  normal  oper- 

i 

ations,  however,  applied  voltages  for  the  sweeps  would  be  obtained  by 

J 

a linear  transformation  on  the  monitor  word  voltages.  In  the  first  1 

satellite,  it  has  been  determined  that  the  nominal  voltages  are  nearly  i 

correct,  and  since  the  actual  voltages  are  noisy  the  nominal  voltages 
are  used  for  the  analysis.  The  electron  and  Ion  nominal  voltages  are 
carried  In  arrays  APEVL  and  APIVL  respectively. 

I ' 

Provision  has  been  made  In  SSIE  for  normal  operation,  that  Is  for  valid 
monitor  word  information.  Over-one  electron  sweep  (10  seconds  in 
length)  there  are  60  monitor  words,  and  over  an  Ion  sweep  there  are  72. 

The  monitor  word  voltages  and  the  corresponding  times  are  collected  In 
the  arrays  EMVOL  and  EMTIM  respectively.  A straight  line  fit  Is  done 
on  this  data  to  obtain  voltage  as  a function  of  time  with  Intercept  A 
and  slope  B.  The  actual  applied  voltages  are  then  obtained  from  the 
times  of  data  words  of  the  electron  and  Ion  sweeps  (In  arrays  XTIME  and 
XTIMI  respectively)  and  the  straight  line  equation.  The  applied  vol- 
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BUFFER  IN  an  1120  word  record 
to  672  CDC  words  of  BF,  L=L+1 


not  ready,  or  EOF 


DO  301  K=l,  224 


Transfer  5 UNIVAC  words  in  3 wds  of  BF,  to 
5 wds,  each  left  justified  with  4 blank 
fill  characters  on  right, of  DEC 


Print  mess! 
than  614  r< 

ige,  "more 
•cords  " 

K=l,  10 

n 

fv_ 

I * (K-l)  *112+1 
M = M + 1 

CALL  WRITMS  (2,  DEC(I),  112,  M) 


CALL  CLOSMS  (2) 


Figure  3.  Program  GWRIT  Flow  Chart 


tages  are  stored  in  arrays  APEVL  and  APIVL  respectively. 


In  SSIE,  lines  43-62,  the  nominal  applied  voltages  for  the  ion  and 
electron  sweeps  are  calculated  and  placed  in  arrays  APIVL  and  APEVL 
respectively. 

Since  each  second  contains  20  data  words,  the  ion  sweep  (duration 
12  seconds)  has  240  data  points  and  the  electron  sweep  (duration  10 
seconds)  has  200  data  points.  Both  sweeps  start  with  the  last  monitor 
word  in  a second.  This  starting  point  then  is  -5V  for  the  ion  sweep, 

and  5V  for  the  electron  sweep.  The  end  of  the  sweep  is  at  12V  in  the 

ion  case  and  at  -8V  in  the  electron.  Thus  in  filling  the  APIVL  and 
APEVL  arrays,  we  take,  in  the  ion  case,  84  points  from  the  240  by 
calculating  voltage  at  the  2nd  point  (the  last  of  the  initial  second), 
then  the  2nd,  5th,  7th,  11th,  14th,  17th,  20th  points  from  the  2nd 

through  12th  seconds,  and  the  2,  5,  8,  11,  14,  17  from  the  last  second 

according  to  the  formula 

APIVL  (KAP ) = 17/240  * (M-l)  - 5. 

Similarly,  we  calculate  the  nominal  electron  voltages  In  APEVL,  using 
the  next  to  last  point  from  the  Initial  second,  the  1,  4,  7,  10,  13, 

16,  17  points  from  the  2-10  seconds,  and  the  1,  4,  7,  10,  13,  16  from 
the  last  second  according  to  the  formula 

APEVL  (KAP)  * -13/209  * (M-l)  + 5. 

The  GWC  1110  computer  has  random  access  mass  storage  files.  These 
files  are  sector  addressable,  each  sector  being  28  words  long,  and  the 
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first  sector  being  sector  0. 


To  mimic  the  GWC  random  access  on  the  AFGL  CDC  6600  computer,  we  use 
the  READMS,  WRITMS,  OPENMS,  and  CLOSEMS,  random  access  mass  storage 
functions.  It  is  possible  to  directly  mimic  the  28  word  addressable 
sectors  in  this  way,  but  since  the  number  of  sectors  at  GWC  is  usually 
large  (on  the  order  of  25000),  and  since  at  AFGL  each  addressable  re- 
cord requires  an  index,  directly  mimicking  the  GWC  RAMS  would  utilize 
large  portions  of  core  for  indeces.  Because  of  this  and  also  because 
all  I/O  at  GWC  is  in  112  wd,  or  4 sector  modules,  the  addressable 
records  at  GL  are  chosen  to  be  112  words  long.  This  would  indicate 
that  all  sector  addresses  of  I/O  blocks  at  GWC  would  be  multiples  of 
4 (-1  since  the  first  sector  is  at  0).  Thus  for  the  record  address 
of  sector  #ISAD,  we  would  divide  by  4 and  add  1.  To  ensure  against 
the  case  where  the  sector  address  is  not  a 4 multiple,  we  read  in  two 
112  wd  records,  the  first  at  the  integer  less  than  or  equal  to  ISAD/4 
plus  1,  and  the  following  one.  Then  we  tabulate  by  the  number  of  sec- 
tors in  the  remainder. 

Line  309  through  line  331  of  SSIE  is  a section  to  calculate  spacecraft 
average  velocity  for  a 1 minute  interval. 

The  ephemeris  group  of  a 1 minute  input  data  set  contains  in  words 
11-13  the  earth  center  coordinate  cosines  for  the  end  of  the  minute, 
and  words  14-16  for  the  beginning  of  the  minute.  The  coordinates  are 
in  UNIVAC  1100  floating  point  form,  and  the  routine  BLKCONV  is  used  to 
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convert  them  to  CDC  floating  point.  FACTA  and  FACTB  are  latitude  cor- 
rected radii  of  the  earth  at  the  beginning  and  end  of  the  minute,  resp- 
ectively. The  distance  travelled  is  approximated  by  a straight  line 
and  the  velocity  is  calculated.  AVEL  contains  average  velocity  for  the 
minute,  and  is  available  for  further  use. 


In  program  SSIE,  the  arrays  AMPEL  and  AMPIO  contain  logs  (base  10)  of 
currents  from  electron  and  ion  sweeps  respectively.  The  data  is  re- 
ceived, however,  in  raw  voltages  between  0 and  5.11  volts  in  digital 
form  from  0-511.  There  is  a linear  relationship  between  raw  voltages 
and  logs  of  the  currents  they  represent.  There  are  two  calibration 
cycles  every  1024  seconds,  which  are  detected  by  the  monitor  word  vol- 
tage level.  During  the  first  calibration  cycle  the  electron  word 

voltages  correspond  to  a current  of  -1  x 10~6  amps  and  ion  word  vol- 
-8 

tages  to  5 x 10  amps.  During  the  second  calibration  cycle  the  vol- 

-9  -11 

tages  correspond  to  -1  x 10  amps  and  5 x 10  amps  respectively. 
Each  calibration  cycle  is  2 seconds  long  and  thus  corresponds  to  14 
sample  points  of  electron  and  ion  data  words. 


In  the  program,  calibrations  are  flagged  by  the  variables  IC1  and  IC2. 
IC1  is  set  at  -1  at  the  start  of  the  program.  When  the  first  cali- 
bration cycle  is  encountered,  IC1  Is  set  to  1.  Otherwise,  IC1  Is  1 
during  the  first  calibration  cycle  and  0 In  between  first  calibration 
cycles.  IC2  Is  1 during  second  calibration  cycles,  and  0 In  between 
second  calibration  cycles.  During  the  calibration  cycles  the  14  elec- 
tron and  ion  voltages  are  collected  in  KC1E,  KC2E,  KC1I,  and  KC2I 
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arrays.  At  the  end  of  the  calibration  cycles,  subroutine  CAL1B  Is 
called  to  average  the  14  voltage  arrays,  and  to  calculate  the  coef- 

t 

ficients  in  the  straight  line  equations 

AMPEL(NEL)  = AINE  + BSLE  * IEL  * .01 
AMPIO(NIO)  = AINI  + BSLI  * ION  * .01 
which  gives  logs  of  electron  and  ion  voltages  respectively,  from  raw 
voltages  IEL  and  ION. 

Several  plotting  programs  have  been  written  to  display  the  ion  and 
electron  sensor  data.  Mode  1 raw  data  plots  display  monitor,  ion  and 
electron  sensor  voltages  versus  time  and  location  information.  The 
program  uses  the  preprocessed  file  directly  as  input.  The  SSIE  program 
has  an  option  to  create  a file  containing  mode  2 (sweep)  data  which  Is 
placed  on  a system  file.  The  M2  plot  program  uses  this  file  to  produce 
plots  of  electron  and  ion  sweeps.  The  Ml  Ion  density  file  creation 
program  LOOK  uses  as  Input  the  preprocessed  file,  and  creates  a file 
containing  approximate  ion  densities  and  ephemeris  information.  This 
file  is  used  as  Input  to  programs  which  plot  ion  density  versus  ephem- 
eris linear  In  universal  time;  Ion  density  versus  ephemeris  linear  in 
magnetic  latitude;  and  a three  dimensional  plot  of  several  orbits  of 
ion  density  versus  magnetic  local  time.  The  plotting  programs  are 
described  in  more  detail  In  Section  4.2. 

Figures  4 and  5 are  flow  charts  of  the  logic  of  program  SSIE.  Figure  4 
covers  the  Input  and  output  of  data  sets,  and  Figure  5 covers  handling 
of  data  within  the  data  sets. 
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Set  up  constant  curves  for  applied  ION  and 
electron  voltages  In  APIVL,  APEVL  arrays 
respectively 


LSW  - 706724B  (Sync  word) 


Read  Input  and  output  directories  to 
core  arrays  DIRI,  GIRO,  respectively 


IPTR  « bits  IP-35,  first  word  of  Input  directory  (pointer 
word  number  of  last  readout)  - add  2 for  Index  In  DIRI 


KTIM  - 0 


Set  Initial  estimates  of  coefficients  for  volts 
Amps  linear  transformation  In  AINE,  BSLE  (electron) 
and  AINI,  BSLI  (Ion) 


I CL  * -1 
NP  - NP+ 


Increment  IPTR,  mod  26 


NE.  27 


Decipher  flag  bits  of  pointer  word  to  IFL 


From  current  pointer  word 

ISAD  ■ bits  4-19  (readout  sector  address) 

IDSTS  ■ bits  20-35  (number  of  data  sets) 


Figure  4a.  Program  SSIE  Flow  Chart 


NBLK  » Number  of  blocks  In  readout 
Check  for  wrap  around 

Read  Information  block  at  ISAO  to  array  INF 
Write  ephemerls  Information  on  TAPE8 


Initialize  Indeces:  II  * 0,  ICAL  ■ 0,  IVOLD  « 0,  KEl  = 0, 
KIO  » 0,  NEL  - 0,  NIO  ' 0,  IOUT  - 0,  I ■ 1,  KEM  - 1 
Set  IBAD,  the  current  block  address. 

Decipher  output  directory  Information. 

Calculate  updated  Information. 


Read  current  block  to  array  INAR 


IWD  » 1569  - 392  * I * wd  address  Index  In  INAR,  1 min  data 


1 

Set  ICDC  * Julian  Date  | 

LT.O.DR.GT.  500  S 

Calculate  average  velocity 


IOUT  « IOUT  + 1 


Set  IWDO:  address  In  NOUTAR,  1 nln  output  data 


I WOO  - 1 


Figure  4b.  Program  SSIE  Flow  Chart 
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Figure  4d.  Program  SSIE  Flow  Chart 
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Figure  5e.  Program  SSIE  Flow  Chart  Detail 
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Figure  5f.  Program  SSIE  Flow  Chart  Detail 
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AMU  « (Nj+N2*16)/(N1+  N2 
HP  = (TE  + TEMP) /AMU 


Figure  5g.  Program  SSIE  Flow  Chart  Detail 
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Subroutine  M2EL  utilizes  mode  2 electron  sweep  data  and  produces  elec- 
tron density  and  temperature,  vehicle  potential  with  respect  to  plasma, 
and  current  at  resting  bias  voltage.  The  electron  sweep  analysis  is 
discussed  in  Reference  1. 

Subroutine  FLEPOW  and  its  associated  subroutines  SEARCH,  FUNCT,  and 
CERF  were  written  by  M.  Smiddy  in  accordance  with  Reference  3.  FLEPOW 
utilizes  mode  2 ion  sweep  data  to  produce  densities  of  the  two  species 
H+  and  0+,  average  ion  temperature,  and  vehicle  potential  with  respect 
to  plasma. 

Subroutine  NTRP  is  for  straight-line  interpolation,  used  by  M2EL. 

FITLIN  is  a straight-line  fitting  subroutine,  also  used  by  M2EL.  The 
CALIB  subroutine  is  called  by  SSIE  to  recalibrate  the  volts — ►amps 
linear  transformation  coefficients. 
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4.2  Plotting  Programs  at  AFGL 

1.  Program  LOOK  prepares  data  files  for  plotting  programs  TXSC,  TXCUT, 
TXCGM.  LOOK  utilizes  the  preprocessed  data  file  as  input,  and  sep- 
arates and  transfers  monitor  word,  electron  word,  and  ion  word  data  to 
tapes  2,  3,  and  4 respectively.  Tape  5 is  written  on  once  per  minute 
with  ephemeris  information.  Figure  6 is  a flow  chart  of  program  LOOK. 

2.  Program  TXSC  plots  monitor  data,  electron  data,  and  ion  data  on  a 
long  strip  chart.  All  three  are  plotted  versus  universal  time  and 
other  positional  information.  The  program  is  flow  charted  in  Figure  7. 
A sample  of  the  plot  is  included  as  Figure  8.  Program  TXCUT  has  essen- 
tially the  same  logic  as  TXSC.  The  difference  is  that  only  Ion  data 

is  plotted,  one  datum  per  second,  and  the  readout  is  compressed  to 
about  a one-foot  length.  Figure  9 is  a sample  compressed  orbit. 


i 


I 
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3.  Program  TXCGM  plots  ion  data  linearity  versus  geomagnetic  latitude. 
The  abscissa  displays  other  positional  data  information  in  addition  to 
geomagnetic  latitude.  Gaps  are  Included  in  the  plot  to  indicate  where 
there  is  no  data  at  the  extreme  poles.  To  achieve  linearity  in  magnetic 
latitude.  It  is  necessary  to  use  ephemeris  information  at  one-minute 
Intervals,  accumulate  differences,  and  Interpolate  within  one-minute 
periods.  Figure  10  Is  a flow  chart  of  the  logic  for  TXCGM,  and  Figure 
11  Is  a sample  plot. 

4.  PERP  Is  a program  to  plot  successive  half  orbits  of  Ion  data  versus 
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Read  card,  start  sector,  number  minutes,  last-1  data  block  processed 
Read  card,  range  of  seconds  to  be  processed 
Read  card,  plot  factor,  type  of  plot 
Initialize  block  number  (NBLK) 

Read  TAPE6,  coordinate  conversion  table 
Set  block  address  I BAD 
Set  KSEOO 

Read  RAF  information  sector 
Set  orbit  information  variables 


Read  RAF,  one  block  of  data  (4  sectors)  to  INAR 


Set  IWD,  index  in  INAR  of  current  minute.  Set  ephemeris  vbls. 


Set  index  (I SEC)  of  current  second  of  data  in  INAR 
Set  time  of  1 second  of  data  LTIM 

For  KSEC  1,  set  zero  time  reference  and  write  header  rec- 
ord on  tape  5. 


Skip  frame  with  LTIM 


Set  plot  time  limit 
Check  against  limit 
Set  geomag.  coordinates 
L-0,  M*0,  N-0 


LT.LSTRT 


Figure  6a.  Program  LOOK  Flow  Chart 


D075 ) J=l,20 


M0D(J-1,3)-1 


L = L + 1 


N = N + 1 


Figure  6b.  Program  LOOK  Flow  Chart 
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Electron  and  Ion  Word  Variation  on  Density  Scale 
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Figure  10b.  Program  TXCGM  Flow  Chart 
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Figure  11.  Compressed  Model  Ion  Density. 

Linear  in  Magnetic  Latitude 


magnetic  latitude  (either  from  -90°  to  90°  or  from  90°  to  -90°)  in  a 
three  dimensional  presentation.  The  objective  is  to  display  qualit- 
ative  changes  in  the  density  "surface"  from  orbit  to  orbit,  over  the 
course  of  about  a day.  Hidden  points  are  not  eliminated,  but  a viewing 

M 

perspective  is  chosen  so  that  few  hidden  points  will  occur.  To  th’s 
end,  a fixed  perspective  has  been  chosen;  but  this  can  be  altered  by 
changing  data  statements  in  the  program.  The  assumptions  about  the 
three  dimensional  coordinate  system  are: 

a)  The  first  semi-orbit  is  in  the  y=0  plane  with  magnetic 
latitude  along  the  x axis  and  densities  along  the  z axis. 
Successive  semi -orbits  are  at  y=l/2",  1",  1-1/2",  etc. 
planes.  The  expected  density  values  are  scaled  between 
z=0,  z=3.  The  semi -orbit  magnetic  latitude  is  scaled 
between  x=0,  x=6. 

b)  The  coordinates  of  the  eye  are  (-10. ",-10." ,10.").  The 
line-of-sight  Is  from  the  eye  to  the  point  (3", 3", 2"). 

The  eye  vertical  is  the  z axis. 

c)  The  projection  plane  is< perpendicular  to  the  line-of-sight 
and  is  5"  from  (3,3,2)  along  the  line-of-sight  toward  the 
eye. 

d)  The  coordinate  system  In  the  projection  plane  is  oriented 
so  that  Its  y axis  is  parallel  to  the  line  from  Z to  F in 
the  3D  system.  The  Intersection  point  has  the  coordinates 
(B,C)  In  the  plane  coordinate  system,  where  Z Is  (3,3,0) 
and  F Is  (3,3,2). 
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Mathematics  of  Transformation 
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The  3D  coordinates  of  the  eye  are  in  array  E.  Those  of  the  sighted 
point  are  in  array  F,  both  arrets  dimensioned  3.  The  distance  from 
E to  F is 

3 1/2 

DE  = £ (E(i ) - F(i)2) 

i=l 


The  unit  vector  along  the  line-of^sight  from  F to  E is 


E(i ) - F ( i ) 
DE 


i = 1 


,3 


So  the  point  along  the  line  5"  from  F is 

P(i ) = F(i)  + 5[(E(i ) - F(i))/DE]  i = 1,...,3 

The  equation  for  the  plane  through  P is 

(x-P(l))  (E(l)-F(l))  + (y-P(2))  (E ( 2 )— F ( 2 ) ) + (z-P(3))  (E(3)-F(3))  = 0 

For  an  arbitrary  point  A(l),  A( 2 ) , A(3),  the  equation  for  the  line  to 
E is  given  by 

x = (E(l)  - A(l))  * T 
y = (E(2)  - A(2) ) * T 
z = (E(3)  ■-  A(3 ) ) * T 

Substituting 

3 

£ [((E{1)  - A(i))T  - P(i ))  (E ( i ) - F(i))]  = 0 
1*1 


r 

i 
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Solving  for  T 


3 

E [P<0  (E(i)  - F(  1 ) )] 

1=1 

T = 

3 

E [(Ed)  - F(1))  (Ed)  - Ad))] 
1=1 


Then  the  coordinates  for  the  intersection  point  are 


X(i)  = (E(1)  - A(1 ) ) * T 1 = 1 3 


Now,  to  select  the  unit  vectors  In  the  plane  (in  3D  coordinates),  we 
solve  for  the  Intersection  point  of  the  line  through  (3,3,0)  and 
(3,3,2)  and  the  projection  plane.  The  equation  for  the  line  is 


x = 3,  y = 3,  z 


Solving  for  z coordinate 


(P(l)  - 3)  (E(l)  - 3)  + (P(2)  - 3)  (E ( 2 ) - 3) 

S(3)  = + P ( 3 ) 

(E ( 3 ) - 2) 


with  S(l)  = 3,  S ( 2 ) = 3 


The  unit  vector  in  the  y direction  in  the  plane  U is 


U(1)  = (P(1)  - S ( 1 ) )/DU  1 = 1 3 


where 


I 
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The  vector  P— >(3,3,2)  is  along  the  1 ine-of-sight,  so  G is  a unit  vec- 
tor perpendicular  to  H where 

3 - P ( 1 ) 3 - P (2 ) 2 - P (3) 

G(  1 ) = , G(2)  , G(3)  = 

DG  DG  DG 

with 

1/2 

DG  = [(3  - P (1 ) )2  + (3  - P(2))?  + (2  - P(3))2  ] 

So  the  x unit  vector  in  the  plane,  V,  is  given  by 
GxU 


where 

V(l)  = G(2)  * U ( 3 ) - U ( 2 ) * G ( 3 ) 

V ( 2 ) = G ( 3 ) * U(l)  - G ( 1 ) * U ( 3 ) 

V ( 3 ) = G ( 1 ) * U ( 2 ) - U(l)  * G ( 2 ) 

The  origin  in  the  plane  system  is 

0(i ) = P(1 ) - B * V(1)  - C * U(i)  i = 1.....3 


where  B = -1,  C = 0.  So  the  x coordinate  of  X is  (X  - <t>)  * V,  or 
3 

W(l)  = C(X(1 ) - 0(1))  * V(i )] 

and  the  y coordinate  is 
3 

W(2)  = L C(X(i)  - 0(f))  * U(i )] 

1=1 
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where  W(l),  W( 2 ) are  the  coordinates  of  X in  the  projection  plane 
coordinate  system. 

In  the  main  program  PERP,  E,  F,  Z are  defined  in  data  statements;  B and 
C are  entered  as  defined  variables;  P,  S,  U,  V and  G are  calculated; 
and  A is  redefined  as  necessary  to  call  subroutine  PLAN.  PLAN  calcul- 
ates T,  X and  W. 

Input  File  Preparation 

Program  TXCGM  is  modified  slightly  to  write  files  called  TAPE7  and 
TAPE9.  TAPE7  contains  all  the  plot  coordinates  for  the  ion  plot  versus 
geomagnetic  latitude  of  one  readout.  TAPE9  contains  for  the  readout 
the  x coordinate  of  the  first  pole  (XFIR),  the  indicator  ( MXMN ) of 
whether  the  pole  is  the  north  (1)  or  south  (-1)  pole,  and  the  total 
number  (KEX)  of  poles  in  the  readout.  TXCGM  so  modified  is  run  succes- 
sively through  the  readouts  which  are  to  be  included  in  the  three  di- 
mensional plot.  After  each  run  of  TXCGM,  an  end  of  file  is  inserted  on 
TAPE7,  but  none  on  TAPE9. 

Program  RECAT  then  uses  TAPE7  and  TAPE9  as  input,  and  produces  TAPE1, 
which  has  plot  information  starting  at  the  first  pole,  and  the  indi- 
cator MXMN  as  the  first  record.  RECAT  is  flow  charted  in  Figure  12. 

Finally  program  INPRP  (Fig.  13)  selects  the  semi-orbits  desired  (card 
data,  1 = first  set,  -1  * set  shifted  by  90°),  and  writes  TAPE6,  which 
contains  the  semi-orbits  to  be  plotted. 


PERP  then  plots  the  semi-orbits.  PERP  is  flow  charted  in  Figure  14. 

The  total  number  of  semi-orbits  plotted  is  read  from  a card.  Figure  15 
is  a sample  three-dimensional  plot. 

5.  PLORB  is  a simple  program  which  directly  plots  sweep  data  as  output 
from  program  SSIE.  A sample  plot  appears  as  Figure  16. 


Figure  15.  3D  Mode  1 Successive  Semi-Orbits 
Ion  Density  vs.  Magnetic  Latitude 
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4.3  SSIE  Output  File  Changes 


f 
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I 
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The  organization  of  the  SSIE  output  file  was  described  in  Reference  2. 

Since  that  report,  there  have  been  several  changes  in  the  structure,  and 

i i 

the  following  will  supersede  the  output  file  description  in  Reference  2.  * 

The  processed  file  (output)  from  the  SSIE  analysis  program  fits  the  < 

description  of  the  DMSP  Special  Sensor  Data  Preprocessed  Circular  File 
Subsystem  Specification  of  AFGWC,  revised  3 July  1976,  with  the  fol- 

» 

lbwing  exceptions.  < M 


1.  The  processed  file  is  22016  sectors  long,  addressable  0-22016. 

Each  data  set  is  one  minute  of  data,  but  is  420  words,  or  15  sectors, 
long.  The  one  additional  sector  is  included  in  the  ephemeris  group 
(header)  making  that  group  60  words  in  length,  rather  than  32  as  in  the 
preprocessed  file. 

2.  As  in  the  preprocessed  file,  one  second  of  data  is  contained  in  six 
36  bit  words,  the  first  being  a time-synch  word,  and  the  remaining  five 
containing  20  nine  bit  words  of  data.  In  the  processed  file,  the 
seconds  themselves  are  in  time  order,  but  one  second  of  data  is  re- 
peated as  it  appears  in  the  input  file*  (Ref.  2,  pp.  73-74). 

* 1 

* The  program  has  the  option  to  include  ion  density  approximation  in 
the  ion  data  words. 


— V 
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3.  The  60  word  data  set  header  (ephemeris  group)  for  the  output  file 
Is  as  follows: 

First  Sector  of  Minute 


Contents 

Latitude,  radians,  O-PI/2,  fp. 
Longitude,  radians,  O-PI/2,  fp. 
Altitude,  NM,  fp. 

Julian  Date,  1-366,  integer 

Seconds  from  0000Z  on  JD,  0-86399,  int. 


First 
data/time 
in  minute 


As  in  words  1-5 


Earth  center  coordinates,  time  of  word  5 


data/time 
in  minute 


Earth  center  coordinates,  time  of  word  10 


Orbit  anomaly  angle,  rad,  fp.  Time  of  word  5 
Orbit  anomaly  angle,  rad,  fp.  Time  of  word  10 
JD  of  latest  sweep  pair  start  time,  integer,  1-366 
Seconds  from  0000Z,  integer  0-86399. 
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tg  designates  time  at  midpoint  of  the  electron  sweep.  If  the  start  of 
an  electron  sweep  occurs  within  a minute,  the  time  In  words  19-20  will 
be  between  the  time  in  4-5  and  9-10. 

21  Electron  density  at  tQ  for  latest  sweep  not  starting  in 

.3 

current  minute,  CM  , fp. 

22  Electron  density  at  tQ  for  sweep  starting  in  current 
minute  if  there  is  one  - otherwise,  same  value  as  wd  21. 
CM"3,  fp. 

23  Ratio  of  electron  density  at  tQ,  to  measured  current  at 
sensor  potential  equal  to  resting  bias  voltage  for  the 
same  sweep  as  in  word  21. 

24  As  in  word  23  for  same  sweep  as  in  word  22. 

25  Electron  temperature  at  tg  for  latest  sweep  - current 
minute  sweep  if  one  starts  in  it,  deg  K,  fp. 

26  Vehicle  potential  with  respect  to  plasma  at  tg,  same 
sweep  as  in  25,  volts,  fp. 

The  numbers  contained  in  the  following  five  words  take  their  values 
at  tg  + 17,  midway  into  the  ion  sweep  following  the  latest  electron 
sweep. 

27  N(l)  H+  density 

28  N ( 2 ) 0+  density 
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Second  Sector  of  Minute 

1 TEMP  Average  ion  temperature 

2 AMU  Average  ion  mass 

3 PH  IS I X Vehicle  pot  wrt  plasma 

The  following.  Scale  Height,  is  computed  at  tQ  + 8.5,  midway  between 
the  two  sweeps. 

4 HP  Scale  Height 

The  following  4 words  contain  coefficients  for  transforming  raw  voltage 
to  applied  voltages  after  a linear  fit  of  raw  voltages  and  times. 

For  electron  sweep: 

APEVL(mv)  = XTIME(mv)  * AVEB  + AVEA 

For  ion  sweep: 

APIVL(mv)  = XTIME(mv)  * AVIB  + AVIA 

Raw  voltages  are  obtained  from  monitor  data  words.  At  present  nominal  vol- 
tages are  being  used,  since  monitor  words  are  noisy.  The  logic  is  included 
for  future  use. 

5 AVEA 

6 AVEB 

7 AVIA 

8 AVIB 


; 


The  following  four  words  contain  coefficients  for  transforming  raw  voltages 
to  logs  of  amps.  They  are  recalculated  every  1024  seconds. 
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For  electron  sweeps: 

AMPEL(NEL)  = AINE  + BSLE  * IEL  * .01 

( 

For  Ion  sweeps: 

AMPIO(NIO)  = AINI  + BSLI  * ION  * .01 

f 

9 AINE  r 

10  BSLE 

11  AINI 

12  BSLI 
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5.  Plasma  Bulk  Motion 

Initiator:  P.  Wlldman 

Project  No:  2311  Problem  No:  4701 

5.1  Introduction 

Disturbances  and  Irregularities  In  the  Ionosphere  and  magnetosphere, 
particularly  In  the  polar  and  auroral  regions,  affect  electromagnetic 
wave  communications,  detection,  tracking,  guidance,  and  early  warning 
defense  systems.  Ambient  plasma  bulk  motion  Is  an  Important  factor  in 
forecasting  Ionospheric  and  and  magnetospherlc  plasma  meteorology  and 
In  understanding  the  processes  leading  to  polar  region  Irregularities 
(Fig.  1).  Such  Irregularities  are  already  known  to  cause  the  scintil- 
lations In  satellite  signals  received  at  ground  stations  (Ref.  1)  and 
In  radar  signals  over  the  horizon  (Ref.  2).  This  work  supports  re- 
search projects  at  AF6L  Involving  In  situ  sensing  by  satellites  of 
plasma  bulk  motion  as  well  as  plasma  densities  and  temperatures.  The 
satellites  concerned  are  S3-2  and  S3-3  flown  at  Ionospheric  and  mag- 
netospherlc altitudes,  respectively.  The  on  board  sensors  are  divided 
Into  three  groups:  the  planar  Ion  sensor  array,  the  non-planar  Ion 
sensor  array,  and  the  spherical  Langmuir-type  electron  sensor. 
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5.2  Data  Processin 


For  data  processing,  three  types  of  data  tapes  are  required,  viz.,  the 
plasma  data  tape,  the  satellite  attitude  (OM)  data  tape,  and  the  satel- 
lite ephemeris  (ORMAG)  data  tape.  The  structures  of  these  data  are 
documented  in  References  3 and  4. 

The  satellite  signals  consist  of  several  patterns  which  have  to  be 
recognized.  They  are  the  telemetry  (TM)  applied  voltage  sweeps,  elec- 
tron sensor  real  sweeps,  electron  calibrates,  ion  calibrates,  and  range 
switchings.  The  flow  chart  of  Figure  2 outlines  the  data  processing 


and  computational  schemes.  Figure  3 displays  some  typical  raw  data  of 
electron  and  ion  currents  measured  on  satellite  S3-3. 


DIRECTION  AND  VELOCITY  * PREFERRED  COORDINATES 

RELATI VE  TO  SPACECRAFT  [MAGNETI C,  ECI  ETC.] 


GMT  25126.16 
RLT  728.27 
LRT  80.45 
LONG  65.64 
M .L  . T .47679  .37 
L.T.  40878.67 
I .LflT.74.60 


25486.16 

307.93 

68.21 

315.09 

20940.21 

14707.96 

74.88 


Figure  3.  Electron  and  Ion  Currents  - Raw  Data 
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environment;  (2)  the  effective  aperture  area  Is  affected  by  the  plasma 
sheath  surrounding  the  sensor  and  the  vehicle;  (3)  there  often  exists 
at  the  same  region  more  than  one  species  of  electrons,  characterized  by 
different  densities  and  temperatures;  (4)  geometry  of  grids  In  the  sen- 
sor; (5)  secondary  electron  emissions  and  also  photo  electron  emissions 
from  the  sensor.  Figure  4 shows  a typical  current-voltage  pattern 
during  a real  sweep  of  the  electron  sensor. 

The  vehicle  potential  variations  In  the  Ionosphere  are  related  to  the 
phenomenon  of  spacecraft  charging  (Ref.  5).  It  has  been  determined 
that  whenever  the  relative  vehicle  potential  Is  greater  than  about  6V, 
the  sweep  lies  outside  the  region  of  validity  of  the  presently  known 
theory,  and  the  data  has  to  be  rejected.  The  plasma  sheath  Is  a re- 
sult of  the  well  known  Debye  shielding  property  of  plasmas  (Ref.  6). 

The  presence  of  multi -electron  species  characterized  by  their  Indivi- 
dual temperatures  and  densities  Is  due  to  electron  flows  In  the  Iono- 
sphere, especially  In  events  such  as  energetic  electron  precipitations 
In  auroral  and  red  arc  regions  and  In  polar  plasma  flow  regions. 
Photoemissions  and  secondary  emissions  are  electrons  knocked  off  the 
sensor's  outer  surface  by  solar  photons  and  energetic  electrons  resp- 
ectively. 
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In  order  to  provide  an  analytically  tractable  model.  It  Is  assumed  that 
only  two  species  of  electrons  are  present;  l.e.,  the  amount  of  a third 
species  Is  negligible.  Mathematically,  the  sensor  current  becomes  ex- 
pressed as  the  sum  of  contributions  of  Ij  and  I2  from  the  two  species 
and  I j from  the  secondary  and  photoemissions. 


REAL  ELECTRON  SWEEP  - VOLTS 


I ( V ) = IjOlj.TpV)  + i2(n2,t2,v)  + I3(V) 


I(V)  = £ 4 r2  EN.e  /— 1 ‘ I-JT 

1=1  1 V 2*™  r* 


r2-R2 


r2  - R2  kT, 


+ exp  (a  + bV) 


where  V = V$  + <t>,  where  <t>  is  the  vehicle  potential,  and  Vs  Is  the  ap- 
plied  voltage  on  the  sensor,  r Is  sheath  radius,  R grid  radius,  £ trans1 
mission  factor,  k Boltzmann  constant,  e electron  charge,  N electron 
density,  T electron  temperature,  and  1 labels  the  species.  Each  sweep 
pattern  is  the  superposition  of  the  sweep  patterns  corresponding  to  two 
single  species.  The  expression  I ( V ) Is  Invariant  under  permutation  of 
the  Indices  1 (=1,2),  and  the  lower  temperature  species  has  a density 
much  higher  than  that  of  the  other  species.  To  determine  the  para- 
meters that  characterize  the  physical  contents  of  the  sweeps,  the 
method  of  steepest  descent  of  Fletcher  and  Powell  (Ref.  7,8)  Is  used 
for  fitting  the  sweep  patterns.  The  parameters  to  be  determined  are 
N^,  Tj,  <p,  r,  a and  b.  Occasionally,  the  results  of  a sweep  have  to  be 
rejected  because  of  over  saturation  of  data,  broken  sweep  due  to  exist- 
ence of  a time  gap  in  which  no  signal  was  received  on  the  ground,  too 
much  noise  In  data,  excessive  vehicle  potential  due  to  spacecraft 
charging,  excessive  photoemission  overshadowing  the  true  signals,  or 
simply  overshooting  in  the  rapid  descent  technique.  Rejection  codes 
are  listed  In  Table  1.  Usually,  the  sweep  results  are  very  useful. 

The  scheme  of  electron  density  correction  computation  at  the  end  of 
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every  real  sweep  Is  shown  In  Figure  5.  Figure  6 displays  the  profile 

of  electron  temperature  determined  by  means  of  S3-3  data.  The  meas-  , 

I 

urement  of  such  a profile  above  1,000  km  has  not  previously  been  mapped 

thoroughly  because  of  the  limit  imposed  by  the  apogee  of  satellites.  , 
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Table  1.  Rejection  Code  used  In  the  computer 
data  processing  of  Real  Sweeps. 


Data  Rejection  Code 

Reason 

1 

Wild  data,  lying  outside  three 
standard  deviations. 

2 

Temperature  density  data  mis- 
matched. 

3 

Temperature  below  600°K. 

4 

Magnitude  of  vehicle  potential 
greater  than  6 volts. 

5 

Steepest  descent  method  does  not 
work. 

Figure  5.  Scheme  of  Real  Sweep  Electron  Density  Computation. 
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Figure  6.  Profile  of  Electron  Temperature  Mapped  by  Satellite  S3-3 
at  Magnetospheric  Altitudes 


5.4  Determination  of  the  Plasma  Flow  Direction 


The  plasma  flow  direction  is  measured  by  means  of  an  array  of  plasma 
ion  sensors  slightly  displaced  from  each  other.  The  technique  is  to 
utilize  the  difference  over  sum  ratio  R..  of  two  measurements  I.  and 

I J I 

I.;  i.e., 

J 


V'j 

R,j  i ♦ i 

i j 


where  the  current  measured  by  the  i-th  ion  sensor  is  given  by  References 
9 and  10. 


I = £ Ne  Aq  [Cos  YQ  Cos  Y$  Cos(PQ  - P$  + Et)  + Sin  YQ  Sin  Y$] 


where: 


= transmission  of  the  aperture, 

= ambient  ion  density, 

= electron  charge, 

= relative  velocity  between  plasma  and  spacecraft, 

= aperture  area, 

= pitch  angle  of  plasma  flow, 

= yaw  angle  of  plasma  flow, 

= pitch  angle  of  ion  sensor, 

= yaw  angle  of  1 on  sensor, 

= angle  subtended  by  the  spacecraft  x*>ax1$  with  the 
projection  of  spacecraft  velocity  V on  the  XZ-plane, 
with  the  satellite  spinning  about  the  y»ax1s. 


..  . 


The  ion  sensors  are  not  symmetrically  mounted  about  the  array  centroid. 
As  a result,  the  use  of  only  two  sensors,  i and  j,  is  not  sufficient 
for  the  determination  of  flow  direction.  An  algorithm  has  been  worked 
out  by  Lai  et  al  (Ref.  11)  utilizing  the  x ratio,  defined  by 

R13  + tan  (2a)  tan  e 

x = 

R13  - tan  (2a)  cot  e 


where  2a  is  the  angle  between  the  normals  of  any  two  adjacent  sensors, 
y and  e,  are  the  skew  angles  about  the  y-axis  and  z-axis  deviating 
from  the  array  center  line,  (see  Table  2).  The  plasma  flow  pitch  and 
yaw  angles  are  given  by  solving  the  following  equations  (Ref.  12) 


1 - x 

R,d  = tan  a tan  ( - Pn  + y - Et) 

i + *R24  Q ' - * 


tan  (Pg  + Et  y) 

| 


tan  a 


[»-»] 

+ xR24] 


(5) 


(6) 


Since  x is  a function  of  R13,  all  four  ion  sensor  measurements  are 
needed  simultaneously  for  the  determination  of  plasma  flow  direction. 
Once  the  flow  direction  is  determined,  the  magnitude  of  flow  velocity 
can  be  derived  from  the  relative  flow  magnitude  measured  by  the  sens** 
ors,  and  the  vehicle  velocity. 
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5.5  Blockage  of  Ion  Flows  and  Non-linear  Function  Fitting 


t 

After  careful  analysis  of  a considerable  number  of  orbits  of  real 
satellite  data,  it  became  clear  that  for  cases  where  the  Mach  number 

i 

of  the  satellite's  motion  relative  to  the  ambient  plasma  is  large  > 

(M  >>  1),  the  current  measured  by  a sensor  is  severely  limited  when  the 
plasma  flow  vector  intercepts  a portion  of  the  spacecraft  surface  reach- 
ing the  surface  aperture.  When  the  flow  vector  does  not  intercept  the 
satellite  surface,  the  current  does  behave  similarly  to  that  detected 
by  an  "unobstructed"  sensor,  as  formulated  by  Lai  et  al  (Refs.  11,12). 


For  the  case  when  the  Mach  number  of  the  satellite  motion  relative  to 
the  ambient  plasma  is  near  or  greater  than  1 (M  Z 1),  the  current  meas- 
ured by  a sensor  is  less  than  that  which  would  be  measured  by  an  "un- 
obstructed" sensor  for  all  directions  of  the  flow  vector. 


To  model  exactly  the  shapes  of  the  obstructions  of  flow  to  the  sensor 
would  require  an  algorithm  that  would  not  be  efficient  for  the  analysis 
of  a large  volume  of  data.  Instead,  the  obstructions  are  approximated 
by  a series  of  wedges  with  wedge  angles  and  <t>2  and  centered  on  the 
sensor's  apertures  (Ref.  13).  As  a result,  the  analytical  expression 
of  the  current  measured  by  a sensor  looking  out  of  the  satellite  spin 
plane  becomes  (Ref.  13): 

oo 

I =1/2  NeAta/ir)1^2  / u i erf  [vS(u  cot  <j>  + q cos  n)3  (7) 

out  J z l z 1 

o 

+ erf  Es/a(uz  cot  4>2  - q cos  £)]j  exp  [-a(  u2  - q cos  e)2]  du2  . 
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where  u,  the  velocity  of  plasma  particles,  is  just  an  integration  vari- 
able here;  q is  the  measured  velocity  of  plasma  flow;  and  9 is  defined 
later. 

Similarly,  the  analytical  expression  for  the  in-plane  sensor  current 
becomes 

oo 

Iin  = 1/2  NeA(a/7r)^2  J „z  | erf  vz  cot  <t>^  + q cos  >?)]  (8) 

o 

+ erf  [v/a(t>z  cot  4>2  - q cos  £)]  } exp  [-a(  v z - q cos  9 )2]  di*z  . 


where 


q cos  d 
q cos  v 

q cos  £ 


(y.  + **)'  • L_, 

(v  + w)  • M, 


2 2 2 
q - (q  cos  v)  - (q  cos  e)  , 


[ v]  = | v 


cosY„  cos  P„ 

q q 

sin  Y_ 


cos  Y sin  P„ 

q q 


r “°v 

cos  Y 

cos 
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sin 
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+ D sin 
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-(cos  Y$  sin  Y$  cos  P$)/x 


-[(cos  Y r cos  P sin  P 1/x 

J j j 


(cos  Y sin  P + sin  Y ) 


JL  f 2 

and  erf (x)  = I exp(-t  ) dt 


These  integrals  can  be  calculated  with  good  accuracy  with  Gaussian 
quadrature.  To  determine  the  plasma  flow  parameters  use  is  made  of 
methods  for  fitting  non-linear  functions  to  data.  Numerical  values  of 
parameters  are  found  by  a search  in  parameter  space  for  the  "solution" 
that  yields  a minimum  in  the  difference  between  the  expected  and  real 
values  of  the  function.  The  search  algorithm  is  based  on  the  method 
of  parameter  optimization  developed  by  Fletcher  and  Powell  (Refs.  7,8). 
Using  the  methods  described  in  this  section,  we  have  computed  the 
plasma  flow  data  for  many  orbits  (Figs.  7,8).  In  particular,  a very 
high  (9.8  km/s)  flow  velocity  of  plasma  was  discovered  to  be  associated 
with  a large  poleward  electric  field  structure  (see  Section  6)  at  the 
plasmapause  during  the  May  1976  magnetic  storm  (Ref.  14). 
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Typical  Ion  Currents  Measured  bv  a Four  Sensor  Array 
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5.6  Remarks  on  Algorithm  Development 

' 

Three  observations  are  made  here  concerning  some  subtle  points  expert 
ienced  in  the  development  of  plasma  flow  algorithms. 

' • > 

The  method  of  steepest  descent  for  the  search  of  roots  is  useful  in 
many  areas  of  science  where  data  analysis  is  used  to  determine  multi- 
ple parameters  in  physical  models.  Instead  of  seeking  rigorous 
solutions  for  the  parameters,  which  may  be  impossible  due  to  the  non- 
linear  nature  of  the  equations,  one  seeks  the  least  square  deviation 
of  the  theoretical  function  from  actual  data.  Starting  from  some 
trial  parameters,  the  computer  algorithm  drives  the  parameters  down 
along  the  steepest  descent  of  the  least  square  function.  However,  due 
to  local  slope  conditions,  overshoots  do  occur.  They  are  normally  not 
a problem,  for  the  parameters  will  eventually  fumble  their  way  back  to 
the  correct  neighborhood.  But  if  the  overshoot  is  such  that  the  para- 
meters have  entered  a region  where  the  equations  are  no  longer  valid 
physically,  the  local  slope  conditions  there  may  mislead  the  root 
search.  To  overcome  this  pitfall,  one  should  build  restriction  cond- 
itions in  the  theoretical  formulation,  or  In  the  algorithm,  so  that 
such  overshoots  into  preestimated  undesirable  regions  do  not  occur. 

- ■ ^ i I | 

The  second  remark  concerns  the  speed  of  descent.  The  function  may 
have  a long  valley  in  parameter  space.  The  algorithm  may  then  happen 
to  shoot  back  and  forth  across  the  valley  numerous  times  during  the. 
descent.  If  one  truncates  the  Iteration  loops  in  an  algorithm  by  im- 
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posing  a certain  number  for  maximum  search  trials,  one  might  end  up 
at  a set  of  parameters  nowhere  near  the  true  solutions.  To  avoid 
this  pitfall,  one  should  scale  the  parameters  so  that  the  long  valley 
deforms  into  a conical  depression.  This  would  greatly  enhance  accur- 

t ■ 

acy  of  solutions  and  minimize  computer  time.  This  is  a potentially 
important  consideration. 

< 

The  third  remark  concerns  the  mathematical  function  erf(x)  as  computed 
by  the  CDC-6600.  When  the  argument  is  negative,  no  diagnostic  indi- 

• t 

cation  is  given  by  the  computer.  It  simply  gives  a zero  value,  so 
that  any  computation  involving  erf(negative)  would  give  a meaningless 

14 

result;  and  the  programmer*-analyst  would  not  know  why.  The  remedy  is 
simply  to  use  erf(-x)  = -erf(x)  in  the  algorithms.  CDC  should  rectify 
this  situation  in  the  future,  or  at  least  warn  users.  , 
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6.1  Introduction 

» 


Electric  fields  are  an  important  aspect  of  the  ionospheric-magneto- 
spheric  region  of  space.  The  plasma  in  this  region,  and  the  associated 
electromagnetic  fields,  form  the  environment  of  the  Air  Force's  satel- 
lite systems,  with  which  they  experience  important  interactions.  The 
reconnection  of  the  magnetic  field  lines  frozen  in  the  solar  wind  ind- 
uces electric  fields  which  are  subsequently  mapped  down  into  the  polar 
magnetosphere  and  ionosphere.  Electric  fields  perturb  the  ambient 
plasma  flow,  affecting  electromagnetic  wave  communications,  especially 
in  the  polar  region.  In  order  to  map  the  global  distribution  of  these 
electric  fields  and  their  associated  currents,  the  Air  Force  scientific 
satellite  S3-2  was  launched  into  a nearly  polar  orbit. 
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6.2  Rotating  Satellite  Boom  Sensor  System 


The  S3-2  satellite  featured  three  pairs  of  deployable,  mutually  ortho- 
gonal booms.  It  had  an  initial  apogee,  perigee.  Inclination,  and  spin 
period  of  1556  km,  240  km,  96.3°,  and  20  seconds  respectively  (Ref.  1). 
The  electric  field  sensors  were  carbon  coated  spheres  of  3.81  cm  radius, 
four  of  which  were  mounted  on  the  ends  of  13.72m  long  triaxial  cable 
wire  booms  deployed  in  the  radial  force  field  of  the  spinning  satellite. 
Two  other  sensors  were  mounted  on  the  ends  of  two  rigid  booms  deployed 
parallel  to  the  spin  axis  with  a tlp-to-tip  11.18m  electrode  separation. 
Data  were  obtained  from  only  one  of  the  wire  boom  pairs  in  the  satellite 
spin  plane  and  from  the  axial  pair  of  sensors.  Experimental  techniques 
for  making  double  probe  electric  field  measurements  have  been  described 
elsewhere  (Refs.  2,3). 

The  frequency  components  of  the  data  from  the  sensors  were  predominant- 
ly at  spin  frequency  and  orbit  frequency.  However,  other  frequencies 
corresponding  to  various  oscillation  modes  may  be  present  during  per- 
iods of  perturbation  to  the  satellites'  motion.  These  effects  are  esp- 
ecially pronounced  during  the  transient  response  periods  following  boom 
deployment,  spin  orientation  maneuver,  etc.  A comprehensive  study  of 
the  mode  dynamics  of  the  rotating  satellite  has  been  fully  documented 
" the  References  4,  5,  and  6,  and  will  not  be  discussed  further  here. 

n the  spin  plane  detected  predominantly  sinusoidal  electric 
<*»  •*.  •*#.  rotational  behavior  of  the  satellite.  The  predoml nant 
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electric  field  measured  is  the  Induced  field  (v.  x B),  where  v is  the 
velocity  of  the  satellite  and  is  the  magnetic  field.  This  induced 
field  is  a relativity  effect  (Ref.  7).  Furthermore,  the  rotating  frame 
of  the  earth  introduces  a Coriolis  effect  (Ref.  8)  (w  x r)  which  also 
has  to  be  subtracted  in  order  to  separate  the  genuine  electric  field, 

» fr°m  the  measured  field  Em  . Thus, 

—true  —me as 


^neas 


—true 


(v 


where  is  the  earth's  rotational  angular  velocity,  and  _r  is  the  pos- 
ition vector  of  the  spacecraft  measured  from  the  center  of  the  Earth. 
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6.3  Data  Structure  and  Processing  Techniques 


For  data  reduction  and  computation,  three  types  of  input  data  are  re- 
quired. They  are:  (1)  electric  field  experiment  data,  (2)  satellite 
sensor  attitude  (OM)  data,  and  (3)  ephemeris  and  magnetic  field  (ORMAG) 
data.  The  structures  of  these  data  are  tabulated  in  the  References  9, 
10,  and  11.  The  electric  field  experiment  data  are  in  the  fifth  file  of 
every  tape  of  the  satellites's  raw  data.  The  experimental  data  contain 
six  channels  of  sensor  signals  with  various  high  and  low  gains.  These 
signals  are  given  in  telemetry  counts,  which  can  be  converted  to  volt- 
ages by  a simple  conversion  factor. 


A central  idea  for  the  determination  of  the  various  geophysical  para- 
meters from  experimental  data  is  the  fitting  of  curves  using  theoretical 
formulations  based  on  physical  considerations.  The  main  program  system 
for  this  project  employs  curve  fitting  techniques  using  satellite  orbit- 
al frequency  u>Q  and  spin  frequency  o>s . Transformations  of  spacecraft  and 
geophysical  coordinates  are  also  crucial  in  the  computation  algorithms. 
Some  of  the  transformation  matrices  are  documented  in  Reference  12. 


6.4  The  DC -Offset 
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Due  to  the  relatively  superior  sensitivity  of  the  instruments  on  this 
satellite,  a new  feature  has  been  detected  in  the  nature  of  the  elec“ 
trie  field  data.  It  is  the  DC-offset,  V(t),  which  is  a very  slowly 
varying  function  with  some  orbital  periodicity.  Correlation  analysis 
using  eq.  (1),  has  been  carried  out  in  order  to  determine  whether  the 
DC-offset  of  the  radial  boom  system  would  mirror  that  of  the  axial 
and  thus  allow  for  rigorous  correction.  I.e., 


V Booms  5-6^ 


a + b 


^ Booms  1-2 


(2) 


This  would  be  the  case,  for  example,  if  the  effect  were  due  to  a mono- 
layer  of  0+  Ions  on  the  spheres.  However,  the  data  does  not  support 
this  relationship.  The  axial  voltage  level  tracks  the  plasma  density 
even  in  a pure  H+  plasma.  This  occurs  even  equatorward  of  the  light 
Ion  trough.  One  view  (Ref.  13)  Is  that  the  offset  V ( t ) is  not  a real 
electric  field  variation  but  some  complicated  function  of  plasma  dens- 
ity and  electric  field.  Thus,  the  offset  V ( t)  Is  an  Instrumental  and 
environmental  effect  and  must  be  subtracted  off  If  true  electric  field 
in  the  ionosphere  is  to  be  determined. 


-true 


^neas 


(v 


x r)  x B - V ( t ) 


- we  - 


(3) 


The  method  of  DC-offset  subtraction  is  as  follows.  For  lower  latitudes 
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( <60° ) , where  the  earth's  magnetic  field  lines  are  closed,  the  ambient 
electric  field  is  vanishingly  small  because  of  the  lack  of  solar  wind 
induced  effects  mapping  down  along  open  magnetic  field  lines.  There- 
fore, any  voltage  that  remains  after  subtraction  of  the  induced  field 
U x B)  can  be  attributed  to  the  sphere's  DC-offset.  By  deleting 
regions  above  60°,  and  by  fitting  a slowly  varying  function  (eq.  4) 
with  orbital  periodicity  u>0*"*  to  this  voltage,  one  can  determine, 
approximately,  the  contact  potential  for  an  entire  orbit. 


V(t)  = £a.  sin  iui  t + £b.  cosi  u>  t + C 
l o 1 o 

i=l  i=l 


In  addition,  there  is  another  instrumental /environmental  effect.  The 
instruments  operate  with  a calibration  sequence  every  512  seconds.  But 
during  and  immediately  after  a calibration  sequence,  the  satellite  data 
show  an  exponentially  recovering  pattern  for  about  40  to  60  seconds. 
These  post-calibration  data  must  be  deleted  before  the  treatment  using 
eq. (4 ) can  be  used. 
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6.5  Program  System  and  Description 


The  main  analysis  is  presently  executed  by  means  of  a highly  integrated 
system  of  programs  as  depicted  in  the  flow  chart  of  Figure  1.  Figure  2 
represents  a more  detailed  flow  chart  of  some  of  the  initial  operations 
shown  in  Figure  1:  calibration,  v x B subtraction,  and  Fourier  or  poly 
nomial  offset  fit.  Similarly,  Figure  3 illustrates,  in  greater  detail, 
some  of  the  final  processing  steps:  main  processing,  geophysical 
transformations,  and  special  techniques  for  fitting  DC  offsets. 

Initial  processing  unpacks  the  raw  data  file,  recognizes  data  words, 
converts  data  to  raw  voltages  by  the  appropriate  amplifier  gain, 
removes  amplifier  bias,  subtracts  the  induction  field,  and  eliminates 
calibration  sequences.  Interactive  processing  displays  the  intermed- 
iate data  on  a CRT  graphics  console,  permits  editing,  and  performs 
various  functional  fits  to  the  edited  data  while  allowing  the  user  to 
vary  pertinent  parameters.  Plots  of  the  final  fit  and  a deck  of  co“ 
efficient  cards  are  output  on  option.  The  program  also  has  a 
Butterworth  filtering  routine  and  an  optional  correlation  analysis 
routine.  A non-interactive  version  can  also  be  used  optionally,  and 
is  presently  the  .adopted  version  in  the  consolidated  automated  program 
system.  Since  one  set  of  sensors  in  the  spin  plane  is  inoperative, 
an  instantaneous  determination  of  Ejt)  is  not  possible.  However,  if 
E(t)  varies  at  a rate  slower  than  the  spin  frequency  of  the  satellite, 
then  the  field  can  be  re solved.  More  specifically,  if  the  frequency 
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spectrum  of  E(t),  defined  by 

2(u>)  = /E(t)e1a)tdt  (5) 

Is  bounded  so  that  |z(w)  | = o forw  then  one  can  determine  E_ 
unambiguously.  For  most  orbits  studied,  this  seems  to  be  the  case. 

The  Main  Program  System  then  determines  £ by  fitting  a function  of  the 
form 

V(t)  = A(t)  sin  (wst)  + B ( t ) cos  (u^t)  + C(t)  (6) 

to  the  spinning  boom  voltage.  A(t)  and  B(t)  are  expressed  in  terms  of 
a predetermined  number  of  Chebyshev  polynomials,  such  that  their  fre- 
quency content  does  not  exceed  Since  the  voltage  measured  Is  given 
by 

V(t)  = E{t)  • R(t), 

where  R(t)  * RQ  ( si n(o>$t ) x + coslu^t)  y),  (7) 

one  sees  that  A(t)  and  B(t)  are  the  x and  y components  of  £,  expressed 
In  a coordinate  system  defined  at  t=0.  The  offset  Is  adequately  given 
by  C(t)  In  (6),  provided  that  It  does  not  vary  rapidly  compared  to 
spin  frequency.  The  third  component  of  the  electric  field  Is  measured 
directly  by  the  axial  boom,  with  offset  subtracted,  as  described  In 
section  6.4. 
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Data 


Figure  1.  Flow  Chart  of  Consolidated  Program  System  for  Electric  Field 
Analysis 
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MAIN  Routine 

• Unpack  raw  data  for  a selected  period 
or  whole  orbit. 

• Call  Subroutine  FILER  in  main  loop. 


Figure  2.  Flow  Chart  for  the  Functional  Description  of  Program 
OBTFILEFIT. 
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Offset  Coefficient 
from  OBTFILEFIT 


Program  EFIELD.  Main 
routine  for  electric 
field  computations. 


Subroutine  HEADER. 
Reads  header  data. 


Subroutine  SELECT4. 
Reads  time  and  data. 


Subroutine  OFFSET. 
Subtracts  offset  from 
axial  boom  data. 


Subroutine  LEASTSQ  and  GFUNC. 

Fits  60  sec.  of  spin  plane  boom 
data  using  the  function 

A(t)  sin  w$t  + B ( t ) cos  w$t  + C 

where  A(t),  B(t)  are  Chebyshev 
polynomials.  Only  the  midsection 
20  sec.  fit  is  used  for  e -field 
computation. 


Return  to  main  routine. 
Coordinate  transformation. 


Generate  output  files  of 
electric  field  and  time. 


Figure  3 • Flow  Chart  for  the  Functional  Description  of  Program  EFIELD. 
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6.6  Presentation  of  Results 

f 

l 

The  electric  field  results  obtained  from  the  EFIELD  program  are  trans- 
formed to  various  geophysical  coordinate  systems,  written  on  a disc 
pack  file,  and  plotted  a strip  charts  and  polar  plots.  In  the  strip 
charts,  the  ionospheric  electric  field  magnitude  and  components  are 
plotted  for  a whole  orbit  in  earth  centered  inertial,  geomagnetic, 
satellite  local  vertical,  or  optionally  solar  magnetic  coordinates. 

A typical  strip  chart  is  shown  in  Figures  4 and  5.  In  the  polar 
plots,  the  geomagnetic  components  of  electric  fields  measured  by  the 
satellite  on  polar  passes  are  presented  (Fig.  6).  A line  is  drawn 
from  each  satellite  position  in  the  direction  of  the  electric  field, 
with  its  length  proportional  to  the  field  strehgth.  The  circles  show 
constant  invariant  latitudes.  The  satellite  track  is  not  a straight 
line  because  the  magnetic  pole  is  displaced  from  the  geographic  pole. 

Magnetic  local  time  advances  counterclockwise  from  midday,  at  the  top. 

The  drift  current  map  (Fig.  7)  is  generated  from  the  electric  field 
polar  plot,  by  computing  the  convective  drift  velocity  from  the  elec- 
tric field  and  a knowledge  of  the  earth's  magnetic  field.  Using  the 
techniques  outlined  in  these  sections,  hundreds  of  orbits  of  satellite 
data  have  been  processed  and  the  electric  fields  mapped.  In  parti- 
cular, the  largest  electric  field  structure  ever  detected  at  low 
altitudes  in  the  ionosphere  has  been  discovered  (Ref.  1).  The  event 
took  place  during  a magnetic  storm  and  the  field  structure  was  in  the 
evening  sector  near  the  ionospheric  projection  of  the  plasmapause. 

The  main  pulse  directed  magnetically  northward.  No  significant  elec- 
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Figure  4.  Strip  Chart  Plot  of  Measurements  from  the  E-Fleld, 

Plasma  Motion,  and  High-Energy  Particle  Flux  Experiments 
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Figure  6.  Geomagnetic  Components  of  Polar  Electric  Field  Mapped  on 
a Background  of  Invariant  Latitude  and  Magnetic  Local  Time 
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6.7  Cross  Polar  Cap  Potential 


For  dawn-dusk  orbits  where  the  electric  field  Is  predominantly  In  the 
forward  direction  of  the  satellite,  the  drift  current  must  be  perpen- 
dicular to  the  trajectory.  For  cross  polar  potential  calculations, 

the  Integral  4>  = / £.ds^  is  digitally  approximated  by  * = £ A<l>4. 

1 1 

More  specifically,  let  Er  be  the  component  of  electric  field  measured 
by  the  radial  detectors  (dipole  booms  1-2)  and  ajs  be  the  distance  which 
the  vehicle  which  has  moved  perpendicular  to  B_  in  time  At 

£ x |v  xB) 

AS  = At 

B2 


where  the  velocity  £ must  be  defined  In  the  same  reference  frame  as  £ 
(i.e.,  rotating  or  not).  Thus,  for  each  half  spin  period  t,  the  pot- 
ential change Ad^is  calculated  by  using  (Ref.  14) 


£ x ( v^  x £) 

A4,i  - “£rj  x 


n 

The  total  <*>  A<t>.  is  then  plotted  at  each  incremental  time  At. 

n 1 i 

Calculated  cross  polar  cap  potential  drops  range  between  20  and  about 
80  kilovolts.  The  implications  of  these  observations  for  our  under- 
standing of  high  latitude  convective  patterns  and  their  relation  to 
the  orientation  of  the  Interplanetary  magnetic  field  have  been  pres- 
ented (Ref.  15).  A typical  <!>  polar  plot  Is  shown  In  Figure  8. 
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6.8  Polar  Electric  Fields  and  Plasma  Convective  Flows  In  High 


Latitudes 


The  phenomenon  of  electric  field  reversals  in  the  high  latitude  regions 
near  the  dawn-dusk  meridian  has  been  known  for  several  years  (Refs.  16, 
17).  The  field  reversals  are  consistent  with  convective  patterns  of 
plasma  flow  in  these  regions.  The  patterns  of  electric  fields  and 
plasma  flows  are  highly  dependent  on  solar  winds  and  therefore  on 
interplanetary  magnetic  fields  (IMF).  Using  satellite  S3-2  results, 
morphological  classification  of  these  patterns  has  been  achieved 
(Refs.  18,19),  and  compared  with  earlier  works  (Ref.  20). 

For  data  analysis,  orbits  were  chosen  from  those  of  the  last  four 
months  of  1976  when  S3-2  was  near  the  dawn-dusk  meridian.  They  were 
selected  if:  (1)  the  satellite  passed  to  invariant  latitudes  > 80°, 
and  (2)  hourly  averaged  values  of  the  interplanetary  magnetic  field 
(IMF)  were  avail aLie.  A total  of  43  northern  and  48  southern  polar 
cap  passes  satisfied  these  criteria. 

For  the  benefit  of  those  not  experienced  in  this  field,  a digest  is 
given  as  follows. 

In  Figure  9,  morphological  patterns  of  polar  plasma  convection  are 
displayed. 
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The  flow  in  the  mid  portion  of  the  pattern  is  anti-sunward  and  that  i 
the  outer  portion  of  the  pattern  sunward,  completing  the  convection. 
The  magnetic  field  £ in  the  polar  region  is  practically  downward  so 
that  the  flow  (E  * B)  is  due  to  an  electric  field  in  the  dawn-dusk 
direction,  with  reversals  occurring  twice  in  the  poleward  auroral 
zones,  as  shown  in  Figure  10. 


Figure  10.  Electric  Field  Reversals  in  Auroral  Regions 
in  a Dawn-Dusk  Cross  Section 
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Thus,  a satellite  would  measure  on  its  dawn-dusk  trajectory  a forward  \ 

component  of  the  electric  field  as  portrayed  in  Figure  11. 
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Figure  11.  Forward  Component  of  Electric  Field  in  a Dawn 
Dusk  Orbit,  as  Deduced  from  the  Pattern  in 
Figure  10 
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This  has  been  a simple  case.  In  general,  various  asymmetries  occur  in 
different  cases  (Ref.  18)  depending  on  IMF  conditions  (Fig.  12). 

Seasonal  effects  have  also  been  observed.  In  particular,  a four-cell 
convection  pattern  (Fig.  12)  with  sunward  flow  in  the  middle  has  also 
been  discovered,  (Ref.  19),  agreeing  with  the  theoretical  conjecture 
of  Maezawa,  (Ref.  21). 
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ELECTRIC  FIELD  COMPONENT  ALONG  SATELITE  .TRAJECTORY  (mV/m) 
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6,9  Anomalous  Resistivity  along  Magnetic  Field  Lines 

Due  to  the  helical  motion  of  charged  particles  along  magnetic  field 
lines  B,  the  probability  of  particle  collision  is  mainly  in  the  dir- 
ections transverse  to  B.  The  particles  move  practically  as  free  part- 
icles along  B and  provide  negligible  resistivity.  This  the  reason  why 
electrical  resistivity  is  substantial  only  in  directions  transverse  to 
B;  and  electrical  fields  parallel  to  B^  are  practically  zero.  However, 
interesting  plasma  interaction  events  occur  in  the  ionosphere  and 
magnetosphere.  On  rare  occasions,  anomalous  resistivity  along 
appears. 

Electric  fields  parallel  to  II  have  been  clearly  found  in  one  of  the 
S3-2  orbits:  orbit  517B  has  been  studied  in  very  careful  detail. 

Large  local  fluctuation  of  magnetic  field  appeared,  due  to  the  pres- 
ence of  strong  current  jz  (i.e.,  BBy/c>x  = j'z)  in  the  local  vertical 
system  (Fig.  14).  Strong  electric  fields,  including  substantial  comp- 
onents along  B,  were  found.  The  event  occurred  after  the  satellite 
entered  the  0+  region,  where  the  energetic  electrons  excited  the  0+  to 
cause  visible  arcing  (Fig.  15).  The  gyrofrequency  of  oxygen  ion  has 
been  identified  in  the  FFT  power  spectrum  of  electric  field  data.  The 
plasma  ion  current  (Fig.  16)  showed  significant  change  in  amplitude 
when  entering  the  0+  region,  because  the  heavy  0+  ion  mass  (compared 
to  the  H+  mass)  renders  a great  difference  between  forward  sensor  ion 
current  and  backward  sensor  Ion  current.  The  particle  energy  counter 
clearly  indicated  an  Inverted  V event  during  this  period.  The  forward 
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ELECTRON  FLUX(crr?-sec-sterf ' PARALLEL  CURRENT  (MA/m2) 


Figure  14.  Disturbance  In  Magnetic  Field  Associated  with 
Parallel  Currents  During  the  Satellite  S3-2 
Traverse  of  an  Inverted  "V"  Event 
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Figure  16.  Electron  density  (A),  Ion  sensor  current  (B),  and  Particle  Energy 
counter  (C)  measurements,  showing  an  inverted  "V"  event  occurring 
shortly  after  the  satellite  enters  0+  region. 


component  of  Z shows  a sharp  rise  as  the  satellite  enters  the  inverted 
V region  and  drops  rapidly  Into  negative  values  as  the  satellite  leaves. 
These  exciting  results  will  be  written  up  as  a paper  Intended  for  pub- 
lication (Ref.  22). 
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